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ABSTRACT 

We quantify the fraction of the cosmic infrared background (GIB) that originates from galaxies 
identified in the UV/optical/near-infrared by stacking 81,250 (~ 35.7arcmin~^) if-selected sources, 
spht according to their rest-frame U ~V ys. V ~ J colors into 72,216 star-forming and 9,034 quiescent 
galaxies, on maps from Spitzer/MWS (24, 70, 160 nm), 77ersc/ieZ/SPIRE (250, 350, 500 M-m), and 
AzTEC (1100 i^m). The fraction of the GIB resolved by our catalog is (67± 16)% at 24 ^.m, (72± 17)% 
at 70 ^m, (76 ± 18)% at 160 ^m, (78 ± 18)% at 250 i^m, (70 ± 15)% at 350 \i.m., (67 ± 13)% at 500 |j,m, 
and (52 ± 9)% at 1100 |i.m. Of that total, about 95% originates from star-forming galaxies, while 
the remaining 5% is from apparently quiescent galaxies. The GIB at A ;$, 200 \im appears to be 

sourced predominantly from galaxies at z ;S 1, while at A <L 200 |a.m the bulk originates from 1 ^ z ;$ 2. 

Galaxies with stellar masses log(Af/MQ) — 9.5-11 are responsible for the majority of the GIB, with 
those in the log(M/M0) = 9.5-10 contributing mostly at A < 250 |^m, and those in the log(M/M0) = 
10.5-11 jointly dominating at A > 350 |J.m. The contribution from galaxies in the log(M/M0) = 9.0- 
9.5 (highest) and log(M/M0) = 11.0-12.0 (lowest) stellar mass bins contribute the least — both of 
order 5%. The luminosities of the galaxies sourcing the GIB shifts from a combination of "normal" 
and luminous infrared galaxies (LIRGs) at A ;$ 160 [a.m, to being dominated by LIRGs at longer 

wavelengths. Stacking analyses were performed using SIMSTACK, a novel algorithm designed to account 
for possible biases in the stacked flux density due to clustering. It is made available to the public at 
www. astro . caltech. edu/^viero/vieroJiomepage/toolbox.html. 

Subject headings: cosmology: observations, submillimeter: galaxies - infrared: galaxies - galaxies: 
evolution - large-scale structure of universe 



^ Email: marco.vieroOcaltech.cdu 

t Herschel is an ESA space observatory with science instruments 
provided by European-led Principal Investigator consortia and writh 
important participation from NASA. 

California Institute of Technology, 1200 E. California Blvd., 
Pasadena, CA 91125 

^ Carnegie Observatories, Pasadena, CA 91101 

3 Hubble Fellow 

** Institute for Astronomy, University of Edinburgh, Royal Ob- 
servatory Blackford Hill, Edinburgh EH9 3HJ, UK 

^ Jet Propulsion Laboratory, 4800 Oak Grove Drive, Pasadena, 
CA 91109 

® NASA Postdoctoral Program Fellow 

■^ Laboratoire AIM-Paris-Saclay, CEA/DSM/Irfu - CNRS - Uni- 
versite Paris Diderot, CE-Saclay, pt courrier 131, F-91191 Gif-sur- 
Yvette, France 

® Institut d'Astrophysique Spatiale (IAS), batiment 121, Univer- 
site Paris-Sud 11 and CNRS (UMR 8617), 91405 Orsay, France 

^ Center for Astrophysics and Space Astronomy 389-UCB, Uni- 
versity of Colorado, Boulder, CO 80309 

^^ Dept. of Physics & Astronomy, University of California, 
Irvine, CA 92697 

^^ Department of Physics, Virginia Tech, Blacksburg, VA 24061 

^'^ Laboratoire d'Astrophysique de Marseille - LAM, Universite 
d'Aix-Marseille &; CNRS, UMR7326, 38 rue F. Joliot-Curie, 13388 
Marseille Cedex 13, France 

^^ Institute of Astronomy, University of Tokyo, 2-21-1 Osawa, 
Mitaka, Tokyo 181-0015, Japan 

^^ UK Astronomy Technology Centre, Royal Observatory, Black- 
ford Hill, Edinburgh EH9 3HJ, UK 

^^ Research Center for the Early Universe, University of Tokyo, 
7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan 

^^ Department of Physics & Astronomy, University of British 
Columbia, 6224 Agricultural Road, Vancouver, BC V6T IZl, 



1. INTRODUCTION 

The cosmic infrared background (GIB), discovered 
in Far Infrared Absolute Spectrophotometer (FIRAS) 
data from t he C o smic Background Explorer ( COBE; 
iPueet et all 119961 : iFixsen et all I1998D . originates from 
thermal re-radiation of UV/optical starlight (and to a 
lesser extent active galactic nuclei, or AGN, emission) 
absorbed by dust grains. The total intensity of this back- 
ground is roughly equal to that of the combined extra- 
galactic UV, optical, and near-infrared backgrounds (the 
"cosmic optical background", or COB) indicating that, 
of all the light ever emitted by st ars, about half ha.s bee n 
absorbed and re-emitted by dust (jHauser fc Dwckll200lt ) . 
While it is thought that the majority of the GIB orig- 
inates from du s ty st a r- forming galaxies (DSFGs; e.g., 
iLe Floc'h et~all l2005t iLagacheet all 120051 : IViero et alj 
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120091 ) ■ how they relate to the sources that make up the 
COB, and what fraction of the CIB is resolvable as op- 
tical sources, is still unclear. 

To definitively answer that question, ideally the CIB 
would be resolved into individual sources and matched 
to optical counterparts, but from the first DSFGs im- 
aged in the submm (e.g., ISmail et al.lll997t iBarger et al.l 
Il998t iHughes et al.lll998l ) it became quickly evident that 
identifying optical counterparts is a non-trivial exercise. 
The angular resolution afforded by single-dish submil- 
limeter observatories results in beams containing multi- 
ple sources, such that in deep observations the spatial 
variation of the sky intensity eventually reachs the so- 
called "confusion limit" (;e.g.. lNguven et al.|[2010t) . This 
situation is made worse by the strong evolution under- 
gone by DSFGs bet ween the present day and z ^ I (e.g., 
iPascale et al.l 120091 ). such that only the bri ghtest '-^ 1% 
of DSF Gs (equivalent to ~ 15% of the CIB; lOliver et al.l 
l2010b[ ) at 250 |J.m is resolvable into point sources. This is 
illustrated in FigurelU where a 0.25 x 0.25 deg^ cutout of 
the SPIRE 250 |J.ni map is overlaid with positions of star- 
forming galaxies with masses between ^ lO^'^'^^^Mo, 
at 2; = 1.0-1.5. The map is smoothed and color-stretched 
to highlight the regions of emission. It is clear that very 
few of the sources are detected individually — the rest 
lie almost exclusively on ridges of faint emission. 

Given ancillary data of sufficient quality, this limita- 
tion can be overcome by stacking. Conceptually, stacking 
is very simple: imagine cutting out hundreds of thumb- 
nails from a map centered on the positions where galax- 
ies are known to be, and averaging those thumbnails to- 
gether until an image of the average galaxy emerges from 
the noise. These positional priors can come in many 
forms, e.g., they could be catalogs of UV, optical, IR or 
radio sources. Note that the output is the average of that 
population in the stacked maps, i.e., there will likely be 
sources whose actual fiuxes are higher or lower. Thus, 
the more homogeneous the sources comprising the input 
list, the more meaningful the stacked flux will be. 

Stacking has been successfully applied to infrared maps 
by numerous groups looking to resolve the CIB with 
resolved sources. Frequently, Multi-band Imaging Pho- 
tometer for Spitzer (MIPS) 24 |a.m sources were used as 
positional priors because bright 24 |^m sources are usu- 
ally correlated with far-infrared/submillimeter emission, 
and because of the large fraction of sources that are re - 
solved in 24 um maps ( ^ 70 %; iPapovich et al.l l2003 ). 
For example. iDole et al.l (j2006| ) showed that much of the 
CIB at 70 and 160 |xm is resolved by sources whose flux 
densities at 24 |j.m are S ^ 60 |i.Jy. Similarly, by stack- 
ing on maps from the Balloon-b orne Large-Apertur e 
SubmiUimete r Teles cope (BLAST: iPascale et al.l l2008D. 



iDevhn et all (I2009D and iMarsden et all ()2009D demon- 
strated that close to the full intensity of the CIB at 
250, 350, and 500 [a.m is resolved by 24 |j,m sources, and 
that roughl y half of the CIB a t 500 |i.m originates at 
z > 1.2. iJauzac et al.l (j2011h estimated the contri- 
bution to the CIB brightness at 70 and 160 |i.m from 
24 |j.m sources by stac king in narrow red shift bins span- 
ning 0^^_^J^05^^ iBerta et all ()201lD . using PACS, 
and iBcthermin et al.l ( 2012b[ ). using SPIRE maps, re- 
constructed number counts to flux densities below the 
confusion limit by stacking with 24 |j,m priors, and found 




Fig. 1.— 0.25° X 0.25° SPIRE 250 \xm cutout of the UDS field, 
which has been smoothed and color-stretched to visually enhance 
the regions of submillimeter emission. Overlaid as white circles 
(diameters of 30") are the positions of star- forming galaxies with 
masses between ~ ]^q9-5— lO.O Mq, in a single redshift slice spanning 
z = 1.0-1.5. Note that very few, if any, of the iC-selected sources 
are resolved in the SPIRE map, but that most lie on ridges of faint 
emission. 

that the integral of their counts resolved 58-74%, and 
55-73% of the CIB at their respective wavelengths. 

Sta cking has been used to address other questions as 
well. IPascale et al.l ()2009l ) , stacking 24 |j,m sources on 
BLAST maps, measured the evolution of the infrared lu- 
minosity density with redshift, finding a significant rise 
in the temperatures and lumino si ties of sources with 
increasing redshift. lOliver et al.l (|2010a[ ) used Spitzer 
"bandmerged" catalogs to measure the mass-dependency 
of specific star formation rates (sSFR) — the star- 
formation rate of the galaxy divide d by its mass - — over 
the redshift range < z < 2. IViero et all ()2012a| ) 
stacked high redshift ma ssive galaxies from the GOODS- 
NICMOS Survey (GNS; IConselice et"aI|[20TTI) on maps 
from PACS at 70 -160 ^m TPoeiitsch et al.ll2010[ ). 870 ^m 
from LABOCA (jWeiB et al.ll2009D . and BLAST, find- 
ing that the bulk of the star formation occurs in disk- 
like galaxies, with a hint that spheroid-like galaxies har- 
bor a low l e yel of star formation as well. Similarly, 
iHilton et al.l ()2012l ) stacked a stellar mass selected sam- 
ple of 1.5 < z < 3 galaxies drawn from the GNS 
on SPIRE maps and found evidence for an increas- 
ing fraction of dust-obscu r ed sta r formation with stellar 
mass. And iHeinis et al.l (l2012f ) stacked ultraviolet se- 
lected galaxies at z ^ 1.5 on SPIRE maps, finding that 
the mean infrared luminosity is correlated to the slope of 
the UV continuum, f3. 

While conceptually simple, in practice proper stacking 
is subtle and not without controversy. For sources that 
are uncorrelated (i.e., not clustered around other galax- 
ies) the technique retur ns an unbiased estimate of the av- 
erage flux density (e.g.. IMarsden et al.ll2009l : IViero et al.l 
l2012a[) . But if sources are clustered — which they in- 
evitably will be — then a bias at some level will be 
present and must be accounted for. 
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Many solutions have been proposed t o addres s this 
problem: some, like iBethermin et al.l (j2012br ) cor- 
rect for boosting with sir nulations. Altern atively, 
IBethermin ct all (|2012bD and lHeinis et al.l (I2012D fit the 
measured correlation function to the excess width of the 
measured sta c ked b eam to estimate a correction. And 
iBourne et al.l ()2011l ) use a median statistic to perform 
their stacking, which is shown to be resistant to biases 
induced by outliers. 

Still another method, dev e loped independe n tly b y 
iKurczvnski fc GawiseI^ J2010[ ). iRoseboom et all ()2010D . 
and lBourne et al](j2012t ). fits for the flux densities of mul- 
tiple (correlated) lists simultaneously, thereby account- 
ing for correlations as a part of the stack. The advan- 
tage of this technique is that it makes few assumptions 
and naturally takes into account the possibility that the 
potential clustering bias may be redshift and luminos- 
ity dependent. Here we build upon this technique to 
simultaneously measure the mean flux densities of galax- 
ies selected by mass and divided into mass and redshift 
bins. 

Our goal is to gain a better understanding of the contri- 
bution to the CIB from galaxies identified in the optical 
and near-infrared. In §|2]we present the data, and in ii l3.1l 
we present our method, demonstrating its effectiveness 
with simulations in § 13.31 We ultimately use it to deter- 
mine the total contribution to the CIB from if-selected 
galaxies (§ 14. 3p . and its dependence on redshift (§ 14. 4p . 
mass and color (^ 14. 5|) . and luminosity (? 14. 6|) . The de- 
pendence of sSFR on these variables will be explored in 
a forthcoming paper (Arumuga m et al. in prep.) . 

When required, we assume a IChabrien (l2003f ) initial 
mass function (IMF) and a flat ACDM cosmology with 
flu = 0.27 4. 17a = 0.726, ffn = 70.5kms-^Mpc-\ and 
0-8 = 0.81 (jKomatsu et al.ll201lD . 

2. DATA 

We perform our analysis o n the UKIRT Infrare d 
Deep Sky Survey (UKIDSS; iLawrence et "all [20?)7h . 
Ultra-Deep Survey (UDS) field, centered at coordinates 
2'^17'"50", -5°6'0". The UDS is the deepest survey un- 
dertaken by UKIDSS, covering 0.8 deg^ in J, H, and K 
to nominal 5ct depths of 26.9, 25.9, 24.9 mag [AB]. Cata- 
logs are based on optical and near-infrared (NIR) data in 
this field, while the maps on which the stacking analyses 
are performed span the mid-infrared to submillimeter. 
Here we briefiy describe the catalog and maps. 

2.1. Optical/Near-Infrared Catalog 

Galaxy positions and redshifts come from a cata- 
log based on the UKIDSS Uhra-Decp Survey (UDS; 
ILawrence et al.ll2007t IWarren et"alll2007l ) Data Release 
8 and supplementary data (Williams et al. 2013, in 
prep.). Source detection, photometry, and spectral en- 
ergy di stribution (SED) fitt ing are described in de- 
tail by IWilliams et al.l (I2009D. with recent u pdates to 
the catalog discussed bv lOuadri et al.l ()2012D . A brief 
summary of the data follows: sources are detected in 
the UDS K-haiid mosai c using Source Extractor v2.5.0 
(jBertin fc Arnoutall996[ ). and fluxes measured in other 
bands from PSF-matched images: u' from archival 
CFHT data BVRi'z' from the Su baru-XMM Deep Sur- 
vey (SXDS; iFurusawa et~ani2008D . and JHK from the 
UDS DR8. The catalog is highly complete to Kab < 
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Fig. 2. — Completeness estimates plotted vs. redshift for star- 
forming and quiescent galaxies in bins of stellar mass. 



24.0 for non-pointlik e obje cts. The de-blend i ng te ch- 
nique of lLabbe et~aI1 (j2006f ) and IWuvts et all (|2007l ) is 
used to extract matched 3.6 |J.m and 4.5 |J.m fluxes from 
deep IRAC imaging in the UDS field. Objects near bad 
pixels in the optical or near-infrared images, or those 
with no optical coverage, are excluded, as well as stars, 
saturated sources, severe blends, or those near the edge 
of the image. These quality checks reduce the catalog 
from 171,392 to 81,250 objects, and the resulting effec- 
tive image area is ~ 0.63 deg^. 

Photometric redshifts and rest-frame colors are de- 
rived by fitting the m ulti-band photometry with EAZY 
(iBramnier et al.ll2008r). S tellar m asses are obtaine d with 
FAST (|Kriek et al.ll2009l) using alChabrici] (l20ql IMF, 
solar metallicity, and iBruzual fc Chariot! ( 2003[ ) stellar 
population models. 

Quiescent and star-forming galaxies are classified 
based on the observed bimodality in a r est-frame color- 
color U — V YS. V — J (hereafter UVJ; IWilliams et al.l 
|2009( ): this technique robustly separates red, dusty star- 
burs ts from red, dust- free, old stellar populations (see 
also IWuvts et al.ll2007l ). In addition to the quality checks 
listed above, sources whose best-fit SED has a high x^ — 
which may be the result of poor photometry, artifacts in 
the images, or of a strong active galactic nucleus (AGN) 
component that is not a part of the EAZY template 
library — are excluded from the parent sample. Approx- 
imately 2,700 sources {^ 3%) exceed this limit; we ex- 
plore the effect that inclusion of this subsample has on 
the total resolved CIB in § 15.31 

Our sample is selected at Kab < 24, where the cata- 
log is essentially 100% complete. We calculate the corre- 
sp onding mass comple teness values in a manner similar 
to iQuadri et al.l ()2012f ). Briefiy, we scale the fluxes and 
masses of galaxies at slightly brighter magnitudes down 
to Kab = 24, and estimate the completeness as a func- 
tion of mass and redshift as the fraction of objects below 
that mass at that redshift. Completeness estimates are 
plotted in Figure [2] 
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TABLE 1 

Nominal and Effective Beam Properties. 



Band 


FWHMnom 


FWHMcff 


Acff 


(urn) 


(arcsec) 


(arcsec) 


(steradians) 


24 


6.0 


6.3 


1.548 X 10-3 


70 


18.0 


19.3 


1.296 X 10-8 


160 


40.0 


41.8 


6.151 X 10-8 


250 


18.1 


17.6 


0.994 X 10-8 


350 


25.2 


23.9 


1.765 X 10-8 


500 


36.6 


35.2 


3.730 X 10-8 


1100 


30.0 


30.0 


3.179 X 10-8 



2.2. Spitzer/MIPS 

We use publicly available Spitzer/MIPS maps at 70, 
and 160 |J.m from the Spitzer Wide- Area Inf r ared Ex- 
tragalactic (SWIRE) survey (jLonsdale et all I2003D in 
the X MM Large-Scale Structure field (xmm-lsS: ISuracd 
|2005[ ): and at 24 |j,m from the Spitzer UDS survey 
(SpUDS; PI: J. Dunlop), DR2. Maps have RMS levels of 
0.5 mJy, 1.8 mJy and 19.9 [iJy, respectively. We note that 
the absolute calibration uncertainties are 4, 7, and 12% 
at 24 , 70, and 160 u.m, respectively (lEneelbracht et ahl 
[2007t iGordon et al.ll2007t IStansberrv et al.ll2007D " 

Following iBcthcrmin et al.l (|2010[ ) , calibration correc- 
tions of 1.0509, 1.10, and 0.98, and aperture corrections 
of 1.19, 1.21, and 1.20, are applied at 24, 70, and 160 ^im. 
Maps are in native units of MJysr^^ (surface bright- 
ness), and are converted to Jybeam-^ for this analysis 
by dividing the maps by the solid angles of the measured 
instrumental point response functions ^ (PRFs), where 
f^boam = f/Io = / PRF dfi/PRFo, and PRFq is the peak 
value. 

Also measured from the PRF is the effective full width 
at half maximum (FWHM) of the best-fit Gaussian, 
which can differ from nominal by as much as '^ 6%. 
This is done by simply finding the 2D Gaussian which 
provides the minimum value when differenced with the 
PRF, within a radius of 1.25 x FWHM (chosen as the 
approximate minimum of the primary lobe). Effective 
area and FWHM for each band are listed in Table [T] 

2.3. Herschel/SPIRE 

We use submillimeter maps at 250, 350, and 500 |i.m 
from the t he Herschel Multi- tiered Extragalactic Survey 
(HerMES: [Oliver et al.ll2012D . HerMES is a guaranteed 
time (GT) key project, and consists of maps of many of 
the well studied extragalactic fields, which are divided 
into tiers of depth a nd area, obse rved with both SPIRE 
(jGriffin et al.l [20Tol ) and PACS (iPoglitsch et al.|[20Tol ). 
The UDS is a level 4 field observed with 20 repeats to a 
depth of 11.2, 9.3, 13.4 mJy (5g) , not in cluding confusion 
noise, which from lNguven et all (120101 ) is 24.0, 27.5, and 
30.5 mJy (5cr) at 250, 350 and 5 00 M.m, respectively. Ab - 
solute calibration is detailed in ISwinvard et all ( 2010( ). 
with calibrati on uncertainties of "^7%. Maps are mad e 
using SMAP (jLevenson et al.ll2010t IViero et al.ll2012bD 
and will be made available via HeDaM'^ ( Roehlly et al.l 
1201 Ih . as part of DR2. 

2.4. AzTEC 

We us e maps at 1 100 M.m observed with th e AzTEC 
camera (IWilson et al.ll2008t iGlenn et allll998[ ) mounted 

^ http : //irsa . ipac . caltech . edu/data/ SPITZER/ docs/mips/ 
^ http://hedam.oamp.fr/HerMES/ 



on ASTE (jEzawa et al.|[200l[200l . The FWHM of the 
AzTEC beam on ASTE is 30 arcsec at 1100 M.m, and the 
field of view of the array is roughly circular with a di- 
ameter of 8 arcmin. The area covered is smaller than 
at the other bands, totaling ~ 0.32 deg^ after cropping 
the noisy outer edge. These data will be presented by 
Ikarashi et al. (in prep.). 

2.5. Color Corrections 

We apply color corrections to convert from the stan- 
dard calibration to the actual measured SED of the 
stacked sources. As the part of the spectrum observed 
depends on the source's redshift, the color correction is 
applied after first finding the best-fit SED in each bin. 
Consequently, each color correction is unique, though the 
difference in any one band across the full redshift range 
is never greater than ^ 10%. The color corrections per 
band, from lowest to highest redshift, are: 0.98-1.07 
(24^m); 0.98-1.06 (70|^m); 0.99-1.05 (160|J,m); 0.99- 
1.08 (250 M-m); 0.98-1.03 (350 ^m); 0.98-1.07 (500 M-m); 
and 0.98-1.07 (1100 ^m). 

3. METHOD FOR UNBIASED STACKING 



As was shown by iMarsden et al.l (j2009t ). stacking is 
formally the covariance between a catalog (or multiple 
catalogs) of positions Ca , containing iV^ sources in pixel 
J, and a map, M. The mean of 7V^ is ^q, which rep- 
resents the average number of sources in catalog a per 
pixel. In the limit that sources are Poisson distributed 
on the scale of the beam, then the covariance is simply 
the mean of the map at positions Cq , so that the average 
flux density of a given catalog is 

^ 3 

where A^pix is the total number of pixels in the map. 

If the catalog (or catalogs) in question is correlated on 
the scale of the beam, /i can simply be replaced with 
the variance, a^. What this does not account for — as 



pointe d out bv e.g.JCharv fc Pop3 (120101) . iSerieant et al.l 
()201Gr ). and lKurczvnski fc Gawiscrl ( 20101 ) — is the possi- 



bility that some other, fainter, and potentially numerous 
sources (or the sources in companion catalogs), may be 
correlated with the sources in that catalog, and that ne- 
glecting them could introduce a bias. 

We now present an algorithm, whose f o rmali sm is 
similar to those of iKurczvnski fc Gawiseil (I2O1O0 and 
iRoseboom et al.l (|2010D . with the difference that only 
samples which could potentially be correlated (i.e., those 
in the same redshift range) are simultaneously fit. In the 
following section we provide the formalism, while step- 
by-step instructions are given in § 13.21 



3.1. Stacking Formalism 

The fol lowing is a generaliza tion of the formalism pre- 
sented in iMarsden et al.l (120091 ) , and is applicable to any 
catalog, including those that are clustered at angular 
scales comparable to that of the beam. For a map, Mj, 
with pixels j, and a set of lists, Sa- 



M-j = Tlj 



a 

Si (^N( - ^i) + . . . + 5„ {Ni - M„) , (2) 
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where the Sa form the complete set of all objects in the 
Universe. 

Note that, unlike in iMarsden et alj ()2009t ). we need 
not assume that N^ be a Poisson-distributed number. 
Furthermore, separate lists can also be correlated, so 
that the covariances between them need not be non-zero. 
However, we still require that the instrumental noise is 
well behaved, i.e., (rij) — 0, so that terms in N^ Uj vanish 
in the sum. 

The amplitudes Sa and N^ in Equation [5] that satisfy 
Mj can be quantified by writing their covariances with 
the map itself: 



Cov(M,^„) = -^^M,iV;^ 



iVn 



Dq 



N^ 



E(^^)'-A^"E^^ 



/V • 



J2NiNi,-fia,Y.^'c. 



(3) 



which can be re-written in matrix form by defining am- 
plitude and covariance vectors: 



S = 



fs,\ 

S2 



\SnJ 



; Cov{M,Na) = 



/CoviM,Ni)\ 
Cov(M, N2) 

\Coy{M, N„)) 



(4) 



and 



'Cov{M,Niy 



N, 



pix 



Cov(M, 7V„) 



\j j I 



/Si' 



\S„ 



(5) 



From the covariances and the n x n matrix which we 
label A, S is then simply 



S = A^i Cov(Af,7V„) 



(6) 



Notice the resemblance that the linear system in Equa- 
tion \5\ bears to that of a least-squares fit 

y = E ""^" = "1^1 + "22:2 + ■ . • + anXn, (7) 

a 

whose residual is given by 



i?^ = E p ~ ("1^1 + "2^ + 



QjjiX^ 



(8) 



In order to minimize this residual, we impose the fol 
lowing set of conditions: 



dai 
da2 



dan 



"^E [y^ ~ ("1^1 + "2a;2 
j 

"^E [y^ ^ ("1^1 + 0,2x1^ 



3 



[V] - (ai^^i 



02X2 



dfrX^ 



^JlXji 



dnXj 



0; 



= 0. 
(9) 



These can be expressed in matrix form as 



3 

3 



33 

E^1^2 E(^2) 
3 3 



3 

Erf' rfJ 

j 



0.2 



/ / -^nUj j I Z_^ ^l^n Z_^ ^2'^n ' ' ' Z^V-^n) j 

\ 3 / \ 3 3 3 I 

(10) 

By comparing Equations [S] and [TUl it is clear that the 

Cq, vector maps into Sa , the yj vector maps into M^ , the 

J2^ayj vector maps into the covariances Cov{Mj,Nl), 

j 

and the x^ vector maps into the N^ (which is mean- 
subtracted). 

Therefore, solving Equation [5] for Sa is equivalent to 
finding the coefficients Oq in Equation [7] via a minimiza- 
tion routine. Specifically, the functional form can be op- 
eratively implemented using known quantities: 

M^Y. ScCa = S,C, +S2C2 + ... + SnCn, (11) 

a 

where we define the C^ as a beam-convolved and mean- 
subtracted version of the N^, . 

3.2. Method in Practice 

Here we present the simultaneous stacking algorithm 
(simstack) used in this analysis, which we also make 
publicly available through an IDL code"^. The simulta- 
neous stack is performed on one map at a time, and one 
group at a time, where groups are defined as catalogs 
which could potentially be correlated. For example, we 
group all lists in the same redshift range together (for 
a total of 8 groups), as we expect galaxies of different 
masses but equal redshifts to be correlated with each 
other. In other words, we assume that galaxies in dif- 
ferent redshift slices are uncorrelated, and can be dealt 
with independently. Then, regardless of the code used, 
the method can be broken into four simple steps: 

Prepare N lists of RA and Dec by group, e.g., we divide 
each group (redshift slice) into 8 lists of mass and 
UVJ color; thus iV == 8. 

Construct N layers, or "hits" maps, one for each list, 
where each pixel in the hits map contains the inte- 
ger number of sources which falls into it. 

www. astro. caltech. edu/^viero/viero_homepage/toolbox.htnil 
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Convolve the N layers with an effective point spread 
function (PSF; we use a Gaussian) whose FWHM 
is equal to that of the the effective instrumental 
beam of the real map'*. 

Regress the N convolved layers with the real map of 
the sky, ideally weighted by the noise^. 

Stacking should be performed on maps in Jybeam"'^. 
Errors can be estimated with a bootstrap technique, as 
described in § 13.41 Systematic errors in the method in- 
clude beam area and calibration uncertainties. Note that 
calibration errors may be correlated between bands of the 
same instrument — an effect that should be accounted 
for when fitting models to stacked flux densities. 

3.3. Testing the Method 

Monte Carlo simulations consisting of 10,000 iterations 
are performed to test the estimator for biases. Two sets 
of simulated maps — one containing Poisson distributed 
(random) sources and the other of realistically correlated 
(clustered) sources — are constructed. Each map is a su- 
perposition of sources of varying mass, with the number 
of sources in each mass bin the same as that of real data. 
The flux densities given to the sources are drawn from 
Gaussian distributions centered on the flux densities of 
the measured mean stacked flux densities, and the width 
of the Gaussian five times that of the uncertainty on the 
stacked values, in order to introduce a significant level of 
stochasticity. 

Clustered simulated maps take their positions from 
the actual positions of the catalog sources, in order to 
properly reproduce higher order correlations. Each map 
is then convolved with a Gaussian kernal approximat- 
ing the instrumental beams, with FWHM values ranging 
from 5 to 35 arcsec. 

We then perform stacking analys es on the maps in two 
ways: 1) with a traditional (e.g., iMarsden et al.l [20091 ) 
estimator; and 2) with the SIMSTACK algorithm. Finally, 
the stacked flux densities are compared to the known in- 
put mean values. The histogram in Figure [3] illustrates 
how the traditional estimator returns an unbiased esti- 
mate of the mean flux (S'stackcd/^'input ~ 1) of an unclus- 
tered simulation, but wit h an uncertainty th at depends 
on the beam size (see also lViero et al.ir2012af ). Similarly, 
Figure m shows the bias vs. source density (in number of 
sources per square arcmin) for clustered simulated maps, 
where the traditional estimator is represented by crosses 
and our simultaneous stacking method by circles. The 
data points are offset vertically for visual clarity. The tra- 
ditional stacking estimator is relatively faithful for beams 
of FWHM < 5", but quickly becomes biased, especially 
in catalogs with many sources. The simultaneous stack- 
ing instead returns an unbiased estimate of the mean flux 
density, with errors a = 1-3%, increasing with increasing 
beam size. 



** If using the actual PRF of the instrument, take care that the 
orientation is correct, and if the field has been viewed at multiple 
angles, that the effective PRF is used. 

^ As dictated by the formalism of § 13.11 take care that the mean 
of the pixels to be fit in each layer equals zero. 



C.8 - 

CO 

^ 0.6 

0) 

N 

D 

E 0.4 

o 
z 

0,2 - 



0.0 




0,8 



0.9 



1.0 



1.2 



Fig. 3. — Test of the traditional stacking estimator on simulated 
maps with randomly distributed (i.e., unclustered) sources. His- 
tograms show the resulting output vs. input flux densities of 10,000 
iterations per beam size, for beams ranging from FWHM = 15-35" , 
and a source density of ~ 2arcmin~-^. The vertical dashed line at 
unity represents an unbiased estimate. For all beam sizes, the es- 
timator is shown to be unbiased, though the errors increase with 
an increased number of sources per beam (i.e., for larger beams). 




Fig. 4. — Test of the traditional and simultaneous stacking es- 
timators on 10,000 clustered simulated maps. Recovered vs. input 
fluxes are measured as a function of source density and beam size. 
The traditional estimator performs well for beams smaller than 
5 arcsec, but quickly becomes biased for bigger beams, particularly 
at higher source densities. The simultaneous stacking algorithm, 
SIMSTACK, on the other hand returns an unbiased estimate in all 



3.4. Estimating Uncertainties 

In addition to measurement errors, we must account 
for potential systematic errors introduced by photomet- 
ric redshifts. To address this potential bias we developed 
an extension of the typical bootstrap technique, hereafter 
referred to as the EBT. Like a typical bootstrap, the 
EBT assembles new bins for stacking from the sources in 
the parent catalog, with the difference that rather than 
simply drawing sources randomly from the original bin, 
all of the sources in the catalog are first perturbed ac- 
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Fig. 5. — Stacked flux densities vs. redshift for star- forming (top row) and quiescent (bottom row) galaxies, in divisions of 
circles represent bins with greater than 50% completeness. Note that the flux densities here have been color-corrected (see § 
and errors are tabulated in Tables |4] and [5] for star-forming and quiescent galaxies, respectively. 
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cording to their redshift uncertainties, and then new bins 
are assembled from the new redshifts and masses. Simu- 
lated redshifts are determined by drawing randomly from 
the redshift probability distribution output by the pho- 
tometric redshift code, EASY. Simulated masses, which 
change depending on redshift of the source, must be es- 
timated as well. However, as estimating a new mass for 
every new redshift for every source would be overly labor 
intensive — particularly considering that most perturba- 
tions from the nominal redshift are rather small — we 
instead use the fact that the mass is a strong function of 
K-h&nd magnitude with some slight additional depen- 
dence on J — K color, and we estimate the new mass 
using the perturbed redshift and the observed magni- 
tude and color. Finally, each simulated catalog is split 
up into bins resembling those of the original stack and 
new stacked flux densities are estimated with SIMSTACK. 
This is done 1,000 times. Measurement and bootstrap er- 
rors are then added in quadrature, though we note that 
the error budget is dominated by the EBT estimates. 

The EBT accounts for the possibility of cross- 
contamination of galaxies across redshift and mass bins, 
in addition to the stochasticity of the catalog members 
measured by the traditional bootstrap technique, thus 
resulting in a more realistic error. We find that the EBT 
increases uncertainties by an average of 22%; where the 
correction in bins with better photometry is less; while 
the correction in bins with poor photometry (i.e., high 
redshift and/or low mass) can be as much as 50%. Note 
that these uncertainties account for both instrumental 
and confusion noise, as well as for any pixel-pixel corre- 
lations that map-making may introduce. 

Lastly, systematic errors from the beam area estimate 
for MIPS maps, as well as calibration uncertainties, must 
be taken into account, particularly when estimating the 
contribution to the CIB from galaxies. These errors are 
estimated empirically by including them into the Monte 
Carlo simulation used to estimate the errors. 

4. RESULTS 

4.1. Stacked Flux Densities 

Stacked flux densities and la uncertainties are shown 
for star-forming and quiescent galaxies in the top and 
bottom panels of Figure \E\ and listed in Tables |4] and [3 
respectively. We find statistically significant signals in 
the majority of the bands and bins, with the noisiest sig- 
nals from bins of lowest masses and highest redshifts, and 



the most robust signals in those of the higher mass bins. 
Also, as expected from the beam size and noise proper- 
ties of the maps, the MIPS 24 |^m and three SPIRE bands 
return stacked flux densities with the highest signal-to- 
noise, while the 70 and 160 |a.m uncertainties are larger. 
The uncertainties at 1100 |J.m are also higher, but that is 
largely a reflection of the area of the AzTEC field, which 
is half that of the other bands. 

We note that the traditional stacking method, as an- 
ticipated from simulations (ii l3.3p . returns systematically 
higher results, with the bias proportional to the strength 
of the clustering, which increases with increasing stellar 
mass. Considering this trend, any method that applies 
one correction for all stacked results should be viewed 
with suspicion. 

Also notable is the significant contribution from galax- 
ies identified as quiescent by their colors, a signal which is 
most prominent from the galaxies in the highest redshift 
bins. Their flux densities in all bands increase steadily 
with increasing redshift, to the point where at z ^ 3, 

they are comparable to those of the most massive star- 
forming galaxies in the sample. It is likely that this is the 
result of misclassification of star-forming galaxies arising 
from low signal-to-noise photometry scattering their col- 
ors into the quiescent plane of the UJV diagram. We 
discuss this and other scenarios in S 15.21 



4.2. Best-Fit SEDs 

Intensities, vl^,, are estimated from stacked flux den- 
sities and plotted in Figure |6] with detections shown as 
circles, while non-detections are shown as 2cr upper lim- 
its. Stacked flux densities trace out the SED of thermally 
emitting warmed dust. While the shape of the SED is 
a superposition of many blackbody e mitters of different 
temperatures (e.g.. lWiebc et all2009t ]. it has been shown 
to be well modeled as a modified blackbody of the form; 



S,^Aj^^B{i^,T), 



(12) 



where B{v,T) is the blackbody spectrum (or Planck 
function) with amplitude A, and (3 is the emissivity in- 
dex. The mid-infrared exponential on the Wien side of 
the spectrum can be replaced with a power-law of the 
form /^ oc J^~", which is added by specifically requir- 
ing that the two functions and their first derivatives be 
equal at the transition frequency. Values for both l3 and 
a in the literature range from 1.5 to 2 (e.g.. iBlain et al.l 
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Fig. 6. — Stacked intensities (vlu) of mass-selected sources in the UDS field divided into bins of mass and redshift. Data at 24, 70, 
and 160 ytra (green points) are from Spitzer /MIPS; those at 250, 350, and 500 \aa arc from Hersc/ie//SPIRE (blue points); and those at 
1100 |jm (red points) are from AzTEC. The error bars represent the la Gaussian uncertainties estimated with the extended bootstrap 
method described in § 13.41 Non-detections are shown as 2a upper limits plotted as downward pointing arrows. The median values of the 
redshift distributions are used to convert flux densities into rest-frame luminosities. The SED is modeled as a modifled blackbody with a 
fixed emissivity index, /3 = 2.0, and a power-law approximation on the Wien side with slope a = 1.9. The solid blue line in each panel 
is the best-fit SED, and the shaded region enclosed by dotted light-blue lines shows the systematic uncertainties due to the width of the 
redshift distribution (interquartile range), which dominates the error budget. 
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TABLE 2 

Pearson correlation matrix for all bands under analysis. 
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0.14 
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0.05 
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0.06 
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0.28 
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0.14 
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0.18 


0.16 
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0.19 
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0.13 
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120031 iDunne at al]l201ll : IViero at al.ll2012aD , while we use 
/3 = 2 and a = 1.9, which represent the mean values of 
the best fits of the individual SEDs. Note that for both 
13 and a, we check that the exact value chosen does not 
significantly bias the result. 

Our SED fitting procedure estimates the amplitude 
and temperature of the above template. For the SPIRE 
p oints, the SEP fitti ng procedure (described in detail 
in lChapin et al.|[2008D takes the width and shape of the 
photometric bands into account, as well as the absolute 
photometric calibration uncertainty in each band. Cor- 
relations due to instrumental noise are estimated and 
accounted for with a Monte Carlo procedure. 

Correlated confusion noise must also be accounted for 
in the fit, as these correlations reduce the significance 
of a detection in single band. That is, not account- 
ing for correlated noise in the measurements of, say, 
the three SPIRE bands, would lead to attributing addi- 
tional weight in the overall fit to these data, potentially 
leading to a b i as. T his is discussed in more detail in 
iMoncelsi et~an ()201l[ ) and IViero et al.l (|2012a[ ). We es- 
timate the Pearson coefficients of the correlation matrix 
for all bands from the beam-convolved maps (Table [2]) . 

Interquartile errors, which reflect the uncertainty in 
dimming due to the width of the redshift bins, are es- 
timated from the distribution of redshifts over the full 
set of simulated catalogs generated as described in § 13.41 
Best-fit SEDs for each mass and redshift bin are shown in 
Figure m as solid blue lines; dotted blue lines and shaded 
regions represent the interquartile errors. 

The best-fit SEDs serve several purposes. The first, 
and our primary purpose, is to infer the contribution 
from galaxies in each bin to the entire CIB (spanning the 
full range of wavelengths of our sample) . Another pur- 
pose is to e stima te infrared luminosities, as described in 
iKennicuttI (I1998D . by integrating the SED between rest- 
frame 8 and 1000 |xm (shown as horizontal yellow dotted 
lines in panels of Figure [5]). These are later used when 
quantifying the contribution to the CIB from galaxies 
classified as "normal" , luminous, and ultra-luminous in- 
frared galaxies (i) l4.6p . Infrared luminosities can be used 
as an indicator of obscured star formation, a topic that 
will be explored in Arumugam et al. (in prep.). Finally, 
best-fit SEDs give a measure of the effective dust tem- 
peratures, which we discuss in § 14.71 For reference, both 
temperatures and luminosities are listed in the panels of 
Figure [51 

4.3. Total Resolved CIB 

We estimate the contribution to the CIB from our K- 
selected galaxy sample (first without correcting for in- 
completeness) by multiplying the emission from each bin 
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Fig. 7. — Top Panel: Contribution to the CIB from 
equally spaced redshift slices. Measured data are shown as cir- 
cles with error bars, while best-fit SEDs are shown at dot-dashed 
lines. The redshift ranges plotted here are 0.0—1.0, 1.0—2.0, 2.0— 
3.0, and 3.0—4.0, are shown respectively in red, brown, green, 
and blue. Also plotted as lavender asterisks is the resolved 
CIB using 24 yun priors from Vieira et al. (in prep.). Bottom 
Panel: Contribution to the CIB in divisions of mass, with star- 
forming galaxies shown as stars and quiescent galaxies (with all 
mass bins combined) shown as circles. Both Panels: Also 
plotted is a selection of measurements o f the total CIB in grey, 
from: Spitzer/Ml PS at 24 yrni (diamond: IDole et al.|[2006l ): IRAS 
at eOvrni (boxes; IMiville-Deschcnes et a .1 120021): g j Piiteer/ MIPS 
at 24, 70 and 160 |^m (asterisks: IBcthcrmin et al.l 120101) : and 
70 and 160 ^mi ( exes: IJauzac et al.l 120111) : Herschel/PACS at 
160 vun (triangle s: IBerta et al.ll2011l) : WHAM at 100, 140, and 
240 |jm (crosses: [Lagachc et al.l 120001 ): and COBg /FIRAS spec - 
tra from ~ 200 to 1200 |am (solid line; ILagache ct~an [2OO0I ). 
Lower limits are from Spitzer /MIPS at 24, 70, and 160 ytm 
llPapovich et al.l 120041: IDole et al.l [20061) : SPIRE at 250, 350, and 
500 um IIBcthermin et al.ll2012bl)-^ and S CUBA at 450 and 850 yxm 
llSmail et al.ll2002l : ISerieant et al.ll2004l) . Our total CIB measure- 
ments and best-fit SEDs are shown in both panels as pink boxes 
and dot-dashed lines, respectively. 

(vl^) by the number of galaxies in that bin, and sum- 
ming them together. Results are tabulated in the second 
column of Table [31 Comparing to measurements of the 
absolute CIB from references listed in the last column 
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TABLE 3 
Total Stacked Intensities 






Band 

(|xm) 


Total Stacking 

(nWm^^sr^i) 


This Worlt 
Completeness Corrected 

(nWm-^sr-i) 


Total Model SEDs 

(nWm-^sr-i) 


Absolute 
Absolute CIB 

(nWm-^sr-i) 


Measurements 

Reference 


24 


1.93 ±0.15 (67 ±5%) 
3.66 ±0.92 (55 ± 13%) 
8.06 ±1.72 (62 ±13%) 
6.92 ±1.45 (66 ±13%) 
4.39 ±1.04 (67 ±15%) 
1.91 ±0.33 (73 ±12%) 
0.07 ± 0.04 (38 ± 23%) 


1.97 ±0.15 (68 ±5%) 
3.68 ±0.94 (55 ± 14%) 
8.23 ± 1.76 (64 ± 13%) 
7.05 ±1.48 (67 ±14%) 
4.51 ± 1.06 (69 ±16%) 
1.98 ± 0.34 (76 ± 12%) 
0.08 ± 0.05 (40 ± 23%) 


1.93 ±0.48 (67 ± 16%) 
4.76 ± 1.19 (72 ± 17%) 
9.76 ±2.42 (76 ±18%) 
8.20 ±1.90 (78 ±18%) 
4.58 ±0.99 (70 ±15%) 
1.75 ±0.36 (67 ±13%) 
0.10 ±0.02 (52 ±9%) 


2.86 ±0.17 
6.60 ±0.70 
12.80 ±6.40 
10.40 ±2.30 
6.50 ±1.60 
2.60 ±0.60 
0.19 ±0.04 


Bothcrmin ct al. (2010) 


70 


Botliermin et al. (2010) 


160 


Berta et al. (2011') 


250 


Laeache et al. (2000~l 


350 


Lagaehe et al. (2000) 


500 


Lagaehe et al. (2000) 


1100 


Lagaehe et al. (2000) 



Note. 



Errors estimated with a extended bootstrap technique described in § 13.41 In parentheses are the percentages of the absolute CIB, as 



measured at 24-160 H-m by MIPS; and at 250-1100 |xm by FIRAS. 

of the same Table, we find that our full sample resolves 
about 60%, 70%, and 40% of the CIB at MIPS, SPIRE, 
and AzTEC bands, respectively. 

Next, corrections for incompleteness are made for sam- 
ples that are more than 50% complete, by dividing each 
bin by its completeness estimate (drawn from Figure [5]) 
before summing. Results are tabulated in the third col- 
umn of Table [S] The choice of 50%, though somewhat 
arbitrary, is chosen because at that point the uncertainty 
in the incompleteness estimate is not yet greater than the 
correction. Note that the uncertainties associated with 
this correction are accounted for as part of the Monte 
Carlo simulation used to estimate the total uncertain- 
ties. We find that the completeness correction adds be- 
tween r^ 1-3% to the total CIB. If we relax the complete- 
ness requirement and correct bins which are as little as 
10% complete, the completeness correction rises to ~ 8- 
15% and the resolved CIB to 97% in the SPIRE bands. 
This suggests that a significant fraction of the remaining 
CIB originates from faint sources. This scenario is dis- 
cussed in more detail in § 15.4.11 In all subsequent CIB 
figures, plotted points are completeness-corrected unless 
otherwise noted, with the total contribution in each band 
plotted as pink squares. 

As previously described, the relationship between 
stacked fluxes from different bands is well approximated 
by a simple thermal dust SED. This can be used to 
roughly estimate the total resolved CIB between bands, 
as well as give us a better handle on the contribution 
from noisy bands (in our case 70 and 160 [a.m). Thus, we 
estimate the total contribution to the CIB from summing 
the best-fit SEDs of Figure IH] (corrected for incomplete- 
ness), and report them in the fourth column of Table [31 

4.4. Contribution to the CIB in Broad Redshift Bins 

Plotted as open circles in the top panel of Figure [7] and 
tabulated in Table [S] are estimates of the contribution to 
the total resolved CIB in four redshift bins. The dot- 
dashed lines connecting the circles represent the equiv- 
alent summed SED fits. Notable is the striking depen- 
dence on the contribution to different bands from dif- 
ferent redshifts — a result of negative if-correction — 
with the z = 0-1 bin dominating the CIB at A ;S 160 |J.m, 

and the z = 1-2 bin the chief contributor at A ^ 250 |J.m. 

This behavior shows the same general t rend as the model 
predictions of iBethermin et all ()201lL top panel of fig- 
ure 11), though it should be noted that the limits of the 
redshift bins are not identical. 



4.5. Contribution to the CIB in Mass Bins 

We estimate the contribution to the CIB in divisions 
of stellar mass by summing completeness-corrected emis- 
sion in rows of Figure [SI with star- forming and quiescent 
galaxies grouped separately, and tabulated in Tables [7| 
and [SI respectively. They are also shown in the bottom 
panel of Figure [7| as color-coded stars and circles (qui- 
escent galaxies are all grouped together), with the dot- 
dashed lines representing the equivalent summed SED 
fits. 

A significant part of the resolved CIB ( ^ 65%) at 

all wavelengths appears to originate from star-forming 
galaxies having stellar masses between log(M/M0) = 
10.0-11.0. The total contribution from more massive 
galaxies is ~ 5%, while that from less massive galaxies 
is 30%. Galaxies from the log(M/M0) = 10.0-10.5 bin 
provide the highest contribution everywhere; while those 
from the log(Af/Mo) — 10.5-11.0 bin provide second 
highest contribution at A ^ 160 |i.m, and those from the 

log(M/M0) = 9.5-10.0 bin provide the second highest 
contribution at shorter wavelengths. Quiescent galaxies 
together, unsurprisingly, contribute very little to the to- 
tal resolved background, of order 5%. 



4.6. Contribution to the CIB as a Function of Calaxy 
Luminosity 

Infrared luminosities, Lfir, are calculated by inte- 
jating the rest -frame SEDs between 8 and 1000 |i.m 
Kennicut"tll 19981 ). In Figure [Sj we plot ipiR as a function 
of redshift, in divisions of mass, with star-forming and 
quiescent galaxies displayed as stars and circles, respec- 
tively. Also shown are 24 |xm stacked results from Vieira 
et al. (in prep.) as lavender asterisks. Infrared galax- 
ies have conventionally been classified by their luminosi- 
ties into luminous (LIRG: L — lO'^^"^^ L0) and ultra- 
luminous (ULIRG: L — 10^^^^^ Lq) infrared galaxy 
classes, illustrated in Figure [Sj as horizontal dotted lines 
and right-handed labels. 

There is strong evolution with redshift in the lumi- 
nosities of all galaxies, though it is notably stronger for 
higher-mass galaxies than lower-mass ones. Particularly 
striking is the evolution of the luminosities of quiescent 
galaxies, since at lower than z ~ 2 they are barely de- 
tectable, while by 2: ~ 3 they are nearly as luminous as 
star-forming galaxies of similar mass. This is likely due 
largely to misclassification of star-forming galaxies aris- 
ing from low signal-to-noise photometry at high redshift 
scattering galaxies into the quiescent part of the UVJ 
plane. Note that although the most massive galaxies at 
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Fig. 8. — Luminosity vs. redshift in divisions of mass, for star- 
forming (stars) and quiescent galaxies (circles). Open symbols 
represent bins with greater than 50% completeness. Results from 
Vieira et al. (in prep.) are plotted as lavender asterisks, which by- 
eye appear to be in rough agreement. Notable is the rapid evolution 
in luminosity of the quiescent population, most likely due to noise 
in the photometry of higher redshift sources scattering star-forming 
galaxies into the quiescent section of the UVJ plane. Despite this 
enhanced luminosity, galaxies identified as quiescent provide only 
about 5% to the total GIB C^l43t. 

high redshift are very luminous, they make up a relatively 
small fraction of the full catalog and thus contribute only 
~ 4% of the total resolved CIB (see § 14.51 and bottom 
panel of Figure [7]). It is possible that some fraction of the 
dust heating is due to active galactic nuclei (AGN), as 
similar behavior has be en seen in individually resolved 
objects at 24 |j,m (e.g.. iMarchesini et al.ll2010D . though 
the optical SEDs of this population on average are not 
indicative of AGN. These scenarios are discussed in ij 15.21 
and §[pl 

In Figure |9] we explore the contribution to the resolved 
CIB from different classes of galaxies: so-called "nor- 
mal", LIRGs, and ULIRGs. Data are constructed by 
summing intensities of bins corresponding to their lumi- 
nosities as determined from the best- fit SEDs (§ I4.2|) . 
Short of 160 |J.m, the contribution from LIRGs and less- 
luminous galaxies is comparable, while at wavelengths 
longer than 160 |^m, LIRGs clearly dominate the resolved 
CIB. The contribution from log(L/LQ) < 11 galaxies 
falls rapidly at wavelengths greater than 160 |J.m, which 
may suggest a diminishing contribution from fainter pop- 
ulations at high redshift — which again is suggestive of 
downsizing — but it could also mean that fainter galaxies 
are simply being missed. The contribution from ULIRGs, 
which, as seen from Figure |8l are located at z ^ 1, peaks 

at longer wavelengths, and is an order of magnitude lower 
than less-luminous galaxies at A ^ 160 |J.m. Note that if 

small numbers of exceptionally luminous sources, ultra 
luminous or hyper luminous infrared galaxies, have un- 
usually high specific star-formation rates, this plot would 
fail to capture their distribution accurately. 
Also overlaid in this figure are predictions from 



Fig. 9. — Contribution to CIB from "normal" galaxies 
(L < 10" Lq), LIRGs (L < IQH-i^ Lq), and ULIRGs (L < 
IO^^^^^Lq). Normal galaxies and LIRGs contribute equally to 
make up most of the intensity at A < 70 \xva, which is more sen- 
sitive to lower redshifts, while at longer wavelengths LIRGs and 
eventually ULIRGs con tribute most to the sign al. Also plotted are 
model predictions from lBcthcrmin et al.l 1)20101 . Figure 13, bottom 
panel), with the LIRG and ULIRG predictions somewhat high. 
Although the model is a simple parametric fit to counts at mul- 
tiple wavelengths, the high estimates for the LIRGs and ULIRGs 
lends weight to the suggestion tha t we a re missing luminous, dust- 
obscured sources in our sample ('^ 15.4.2^ . 



iBethermin et al.l (|2011L bottom panel of Figure 11), a 
parametric backwards-evolution model fit to counts at 
multiple wavelengths. The general trends are well repro- 
duced, while in detail, ULIRGs and to a lesser extent 
LIRGs fall short of model predictions. As we discuss 
§ 15.4.21 this suggests that highly d ust-obscured gala xies 
may be missed by our catalog (e.g.. lDev et al.|[r999() . 

4.7. Average Temperature Evolution for Star-Forming 

Galaxies 

Temperatures of the best-fit SEDs for star-forming 
galaxies are shown in Figure 1101 There appears to be 
a mild decrease of temperature with increasing stellar 
mass (i. e., increasing luminosity) — which was also seen 
by e.g., lElbaz et al.l ()201GD — but otherwise the tem- 
perature dependence is dominated by the redshift of the 
bin. 

Also p lotted is a selection o f measurements from the lit- 
erature. iPascale et al.l (l2009l pink exes) is based on best- 
fit SEDs to stacked flux densities — much like ours — 
but with 24 |xm priors on BLAST maps, and are mostly 
consistent with our flndings. The rest of the data rep- 
resent averages in redshift bins of resolved star-forming 
galaxies. Though consistent with each other, this group 
of data is systematically higher than our temperatures 
at all redshifts. This is not at all surprising as indi- 
vidually detected submillimeter sources, particularly at 
higher redshifts, tend to be exceptionally luminous. 



12 



Viero, M. P. et al. 



45 



CD 40 

D 
O 
(U 

a. 35 
E 

I— 

I 30 



"S 25 

(U 



- 


1 1 1 1 


■6-ThIs study: 109(1^15^/^5,) = 
■ft' This study: log(MsF/Mo) = 
^This study: log{MsF/MQ) = 
^Thls study: log{MsF/MQ) = 
X Pascale et al. (2009) 
XAmblord et ol. (2010) 
+ Elboz et ol. (2010) 
Symeonidis et ol. (2013) 


11.0-12.0 - 
10.5-11.0 . 
10.0-10.5 - 
9.5-10.0 . 




a' 






I '■ 






'- ■ - 


, >: 1 I 


A ■ 


J 






i - 






■— ^ - 




)f X 


3 


V- 


w 


A 




"1% 


/^^ 


Ff 1 ^ 


- 


:J 


ySi 


.■■ 


■4 






_ 




■ 1 1 1 i . 


; 



20 



Redshift 

Fig. 10. — Average temperature of the best-fit SED vs. redshift 
for star-forming galaxies in bins of stellar mass. Open symbols 
represent bins with higher than 50% completeness. Though there 
appears to be a mild reduction of temperature with increasing stel- 
lar mass at any given redshift, the temperatures are predominantly 
a function of th e redshift of the sou rces. Also pl otted arc mea- 
surem ents from IPascalc ct al.l (120091. pink exes) ; lAmblard ct al.l 
II2OIOI. purple asterisks ); lElbaz et al.l II2010I . brown crosses); and 
ISvmeonidis ct ahl I I2013I . yellow triangles). 

4.8. Redshift Distribution of the Resolved CIB 

The redshift distribution of the resolved CIB emission, 
d{vl^)/dz, is measured by summing the completeness- 
corrected intensities of all mass bins separately in each 
band, and is plotted in Figure [TT] The peak intensities 
shift from z ~ 0.5-2 with increasing wavelength, indica- 
tive o f the peak of star format i on occurring atz ~ 1-2 
(e.g., iHopkins fc BeacomI [2006t IWang et al.ll2012D . and 
the redshift sensitivity of the different bands. 

Also plotted are predictions for the entire CIB 
from a selection of rec ently published models. The 
iBcthermin ct al.l ()2012aD model in particular, which is 
bas ed in part on stacking m easurements with 24 |xm pri- 
ors (jBethermin et al.]l2012b[ ). describes the measured dis- 
tribution extremely well. However, we caution that since 
the completeness is a strong function of redshift (see Fig- 
ure [2]), it is not expected that the remaining CIB be dis- 
tributed evenly in redshift. Still, this measurement has 
multiple applications, including providing constraining 
power for halo models of future CIB anisotropy mea- 
surements. 



5. DISCUSSION 



5.1. 



Incompleteness Bias 

The SIMSTACK technique returns an unbiased estimate 
of the flux density when all potentially correlated sources 
are accounted for. When sources are missing, either be- 
cause they are too dust-obscured or intrinsically faint to 
detect, their unaccounted-for flux may potentially bias 
the result. While the UDS is an exceptionally deep sur- 
vey, with greater than 90% completeness at z ~ 3.5 
and log(M/M0) = 10.5, it is possible that the absence 
of lower-mass galaxies, particularly at higher redshifts. 
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Fig. 11. — The redshift distribution in different bands of the re- 
solved CIB emission, d(vli,)/ dz, in roughly-spaced redshift slices. 
Over plotted are mode ls from FV^liantc et al.l (12 0091. dashed g ray); 
IBethermin et al.l (120111. solid yel l ow); IBcthermin ct al.l (I2012ai . dot- 
dashed blue); and IVicro cTal] l l2012bl . solid brown at 25 0, 350 
and 500 |Jm); as well as published data from I Ja uzac et aT 
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2012, 



i) 



lavender crosses at 70 and 160 |am). The IBcthermin et al.l 
model, which is determined from stacking 24 ytm priors, fits remark- 
ably well, suggesting a strong correlation between our different sets 
of catalogs. 

could affect the results. We explore the severity of a 
potential bias with the following simulation. 

Restricting ourselves to star-forming galaxies, we first 
simultaneously stack all mass bins at a single redshift and 
record their fluxes. Next, we remove the lowest mass bin 
and stack again, comparing the fluxes in the remaining 
bins to those of the initial stacked fluxes. We repeat this 
until only the highest mass bin remains. Results are just 
as we might expect: the removal of bins increases the 
stacked fluxes of the remaining bins in a systematic way, 
from negligibly little for one missing bin, to as much as 
3-22% for the highest-mass bin alone; and that the level 
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of the bias increases with beam size. For the case of two 
mass bins removed — roughly equivalent to the scenario 
of our highest redshift bin — the bias ranges over 0.2- 
4.0%, from smallest to largest beam size. 

What docs this mean for resolving the background? 
Although missing sources boost the fluxes of those re- 
maining, the sum of their missing flux counteracts the 
boosting, resulting in a marginal bias. Our simulations 
confirm this, with the removal of up to three bins having 
negligible impact. 

5.2. Contribution from Quiescent Galaxies 

Here we briefly discuss the small contribution (~ 5%) 
to the total resolved GIB from galaxies classified as qui- 
escent. The majority of this contribution originates 
at z ^ 2, in contrast with the rest of the GIB, which 

comes mostly from star forming galaxies at z = 0-2. 
Indeed, at z ^ 2 the most massive quiescent galaxies 

(M = 10""^^ Mq) are ULIRGs, and the next most mas- 
sive (M = 101°-" Mq) are LIRGs. 

The most obvious explanation is simply the misclas- 
sification of star-forming galaxies as quiescent galaxies. 
Although our ability to classify galaxies by their colors 
is robust at low redshifts, at z ^ 2.5 the contamination 

due to measurement errors becomes significant. 

On the other hand, measuring stacked infrared lumi- 
nosities as high as these must require that either most of 
the galaxies at these redshifts are misclassified, or that 
the subset of them that are misclassified must be ex- 
tremely luminous. Aside from measurement errors, it is 
worth noting that there are possible physical explana- 
tions for the FIR emission. Because the star formation 
in quiescent galaxies can only have shut off very recently 
at these redshifts, there may still be a strong interstel- 
lar radiation field heating the dust. The galaxies may 
also have some heavily-obscured star-formation activity, 
although the level of activity would have to be low — 
or the level of obscuration high — so that the contribu- 
tion of young stars to the rest-frame J-band luminosity 
is small, in order for the galaxies to still be classified as 
quiescent (see ij2.ip . Another possibility is that there 
is FIR emission associated with AGN activity, although 
dust heated by AGN would be expected to have high 
temperatures and should lead to increased 24 |xm fiux, 
which we do not observe. 

It may well be that some of these supposedly quiescent 
galaxies are actually infrared l uminous; a s c enario which 
is supported by the findings of IVicro ct al.l ([20123), who 
found hints of submillimeter flux for a sample of massive 
early- type galaxies at z ~ 2; though their sample is clas- 
sified by structural parameters rather than by colors as 
done here. 

5.3. High x^ o-nd Active Galactic Nuclei 
Of the roughly 2,700 sources with poor SED fits (de- 



scribed in § 12. 1|) , approximately half of the fits are signif- 
icant hHrngroved with the addition of an AGN template 
(e.g.. lAssef et al.ll2010[ ). We explore this class of sources 
by stacking them separately from — but simultaneously 
with — the star-forming and quiescent galaxies. We find 
further evidence for AGN, since: i) the SEDs are on aver- 
age '^ 10 K hotter; and ii) the 24 |xm flux density exceeds 



that of the other classes of similar mass. The contri- 
bution to the GIB from this class of sources is ;S 4%, 

and is distributed over the same redshift range as the 
star-forming galaxies. 

While intriguing, fully understanding the contribution 
from AGN is complicated, as the threshold chosen to 
classify a galaxy as an AGN is somewhat arbitrary (e.g., 
lAssef et al.ll201Gr ). It is worth noting that because AGN 
emission should exist at varying levels in all galaxies, 
the total contribution to the GIB attributable to AGN 
may in fact be significant. A full treatment of this class 
of sources requires exploring this threshold with better 
fits using a full set of AGN templates, and is beyond 
the scope of this paper, but will be explored in detail in 
Moncelsi et al. (in prep.). 

5.4. Where is the Missing GIB? 

Based on our best-fit SEDs, we resolve between (52 ± 
9)% and (78 ±18)% of the GIB, with the lowest fractions 
resolved at the three shortest (MIPS) and the longest 
(AzTEG) wavele ngths. This is simila r to the resolved 
fraction found bv lMarsden et all (j2009f ) — although they 
likely suffered from a bias due to neglecting cluster- 
ing within the large BLAST beams (jPascale et al.ll2008D 
— and by Vieira et al. (in prep.) and iBcthcrmin et afl 
()2012b[ ) at far-infrared/submillimeter wavelengths. 

Between 50 and 20% of the GIB thus appears to be 
unaccounted for by our if -selected sample. It is notable 
that the measurements for the absolute GIB are them- 
selves uncertain at the ~ 25% level, so that formally our 
measurements could be said to be consistent with the 
GIB in any one band. However, considering that they 
fall below the nominal GIB in every band, the quoted 
percentage should be a fair estimate. 

Besides the possibility that the absolute GIB is not en- 
tirely of extragalactic origin, the two most obvious can- 
didates for the missing fraction are: 1) a large number 
of low-mass, intrinsically faint sources undetected by the 
UDS survey; and 2) IR luminous, potentially massive, 
but unusually dust-obscured sources also missed by the 
UDS selection. We will now discuss these possibilities. 

5.4.1. Low-Mass Faint Sources 

Our sample is selected at K^^ < 24, reaching a source 
density on the sky of ^^ 36arcmin-^. Fainter, lower- 
mass sources missed by the catalog certainly exist, but 
could they contribute enough intensity to make up the 
remaining GIB? 

If we assume that there are as many undetected sources 
as there are detected, and that at, say, 500 |J.m, their 
fiux densities are 0.1 mJy at all redshifts (not unreason- 
able given negative iiT-correction, e.g., see Figure[5|), that 
would add '^ 0.25nWm^^ sr^^, i.e., another '^ 10%, to 
the GIB. 

Also note that in § 14.31 we demonstrated how aggres- 
sively correcting for completeness brought the GIB to 
within 5% of the nominal total, but were held back due 
to large uncertainties in the completeness estimate at 
that level (see discussion in § 12. ip . Thus it is clear that a 
significant fraction — perhaps more than half — of the 
missing GIB is attributable to missing faint sources. 



14 



Viero, M. P. et al. 



5.4.2. Dust Obscured Sources 

Attenuation of UV and optical light in the dustiest 
sources may make them difficult to detect at those bands, 
though of course in the infrared t h ey ca n be quite lu- 
minous. For example, iDev et ahl (jl999| ). observing a 
log(L/LQ) = 12.8 ULIRG at z = 1.44 with the Wide 
Field Planetary Camera 2 (WFPC2; iHohzman et all 
Il995[ ) on the Hubble Space Telescope with the F814W fil- 
ter (A = 7930 A) measured a magnitude of 24.6±0.1 [AB]; 
i.e., this source would probably not have made our cut. 
A hint that obscured sources might be going undetected 
comes from Figure [9l where the res olved fraction of 
ULIR Gs is in disagreement with the iBethermin et al.l 
()201lt ) model, while less-luminous sources appear to be 
in quite good agreement. Since the sum of the three col- 
ored curves resolves the full CIB, it may well be that this 
missing component is rather significant. 

Likely better tracers of highly dust-obscured star- 
forming galaxies are 24 |j,m sources, which, as already 
mentioned, are used by several groups as positional pri- 
ors. Vieira et al. (in prep.) find that their sample resolves 
2.0±0.2nWm~^ sr~^, with which we agree very well. A 
considerable fraction of their emission is shown to origi- 
nate from galaxies hosting AGN. As we discussed in iJ IS 31 
the inclusion of high x^ sources increased the stacked 
24 i^m flux density more than at other bands. However, 
although we resolve about the same amount of the CIB, 
it does not mean we are resolving the same sources. It is 
possible that other faint AGN are being missed, partic- 
ularly at high redshift, where the AGN contribution in 
Vieira et al. (in prep.) dominates. 

Lastly, it is likely that some contribution from "star- 
burst" galaxies — irregular galaxies with very high sSFR 
compared to the main sequence — may missed by the 
mass selection. The potential thus exists that we are 
missing a significant contribution from dust-obscured 
galaxies. In the future, identifying 24 p.m and other 
sources missed by our if-selected catalog — and stacking 
them simultaneously with our optically selected sources 
— would be a natural extension of this work. 

6. SUMMARY 

A complete understanding of galaxies and galaxy evo- 
lution requires characterizing their full SEDs, a goal 
which has to-date been technically challenging, as galax- 
ies in the submillimeter are mostly unresolved. And for 
those that are resolved, large instrumental beams makes 
counterpart identification laborious. Here we attempted 
to statistically connect galaxies selected in the optical K- 
band to those that make up the CIB by stacking sources 



grouped into bins of mass, redshift, and color — present- 
ing and making public our own SIMSTACK algorithm. 

We found that between (52 ± 9)% and (78 ± 18)% of 
the total CIB is resolved by our catalog, and that the 
bulk originates from z = 0-2 and log(M/Mo) = 10.0- 
11.0. Sources at higher redshifts contribute more at 
longer wavelengths, a consequence of the peak of the 
SED redshifting into the submillimeter/millimeter bands 
(i.e., negative X-correction). Higher-mass LIRGs and 
ULIRGs are seen to contribute at longer wavelengths 
(i.e., higher redshift), which shifts to lower mass LIRGs 
and "normal" star-forming galaxies at shorter wave- 
lengths, an indication of downsizing. 

Galaxies identified as quiescent by their colors have 
very little emission at low redshift, but show a rapid evo- 
lution beyond z > 2 most likely indicates a limitation in 
the UVJ color selection at higher redshift due to noisy 
photometry. Lastly, galaxies whose photometry is poorly 
fit by standard SEDs appear to have a clear association 
with AGN; having comparable luminosities but higher 
temperatures as star-forming galaxies of the same stellar 
mass. 

Our work pushes the boundaries of characterizing 
low mass and high redshift optical galaxies in the far- 
infrared/submillimeter, but nevertheless, more can be 
done. Though the UDS is quite deep, deeper and wider 
fields will soon exist. That, combined with future ancil- 
lary data over large areas and to longer wavelengths, and 
with deeper near-infrared imaging and 24 \xni catalogs, 
should allow us to reach the coveted goal of resolving the 
full cosmic background. 
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TABLE 4 
Stacked Flux Densities of Star-Forming Galaxies 



(urn) log(i\I/M(n) = 9.0-10-0 log(M/M(7)) = 10.0-11.0 



z = 0.0-0.5 

log(M/MQ) = 11.0-12-0 



log(M/M(7)) = 9.0-9.5 loe{M/M(^) = 9.5-10.0 



24 


(6.9 ± 0.7) X 10-^ 


(2.2 ± 0.2) X 10~^ 


(5.0 ± 0.3) X 10- 


70 


(7.6 ± 1.3) X 10~1 


(2.6 ± 0.2) X lo" 


(5.2 ± 0.5) X lo" 


160 


(2.2 ± 0.6) X lo" 


(S.l ± 1.0) X lo" 


(1.5 ± 0.2) X IqI 


250 


(1.8 ± 0.2) X 10*^ 


(5.5 ± 0.4) X lo" 


(1.3 ±0.1) X 10^ 


350 


(1.1 ± 0.2) X 10*^ 


(3.1 ± 0.3) X lo" 


(7.3 ± 0.6) X 10° 


500 


(4.6 ± 1.9) X 10"^ 


(1.6 ± 0.3) X lo" 


(3.6 ± 0.5) X lo" 


1100 


— 


(1.0 ± 0.5) X 10^1 


(1.7 ± 0.7) X 10- 



(4.0 ± 0.5) X 10-^ 
(4.2 ± 0.8) X lo" 
(1.4 ± 0.3) X 10^ 
(1.5 ± 0.1) X 10^ 



(9.4 ± 1.1) X 10 
(4.3 ± 1.0) X 10' 
(4.8 ± 1.5) X 10" 



(7.0 ± 0.0) X 10^^ 
(1.6 ± 28.0) X 10^ 
(2.8 ± 26.0) X 10^ 
(3.5 ± 0.2) X 10^ 
(1.5 ± 0.2) X 10^ 
(2.1 ± 1.7) X lo'^ 
(3.6 ± 0.7) X 10-^ 



(urn) log(Af/MQ) = 9.0-10-0 log(-M/M0) = 10.0-11-0 



z = 5-10 

log(Af/MQ) = ] 



0-12-0 log(Af/M0) = 9-0-9-5 log(A//M0) = 9-5-10-0 



jal 



24 

70 
160 
250 
350 
500 
1100 



(1-5 ± 0-1) X 10" 



(5-3 ± 0-8) X 10" 
(4-1 ± 0-9) X 10" 
(2-5 ± 0-8) X 10" 



(8.8 ± 0.4) X 10-2 


(2.1 ± 0.1) X 10- 


(7.4 ± 0.7) X 10-^ 


(1.2 ±0.1) X 10 


(1.9 ± 0.4) X lo'^ 


(6.8 ± 0.7) X 10 


(2.7 ± 0.1) X lo" 


(6.7 ± 0.3) X 10 


(1.8 ± 0.1) X lo'^ 


(4.8 ± 0.3) X 10 


(6.7 ± 1.2) X 10-^ 


(2.1 ± 0.2) X 10 


(1.4 ± 2.0) X 10-2 


(1.9 ± 0.3) X 10- 



(3-8 ± 0-1) X 10-1 
(2-2 ± 0-2) X 10*^ 
(1-1 ± 0-1) X 10^ 
(1-1 ± 0-1) X 10^ 
(8-6 ± 0-5) X lo" 
(4-8 ± 0-4) X lo" 

(3-2 ± 0-7) X 10-1- 



(2-2 ± 0-3) X 10-1 
(1-1 ± 0-5) X lo" 
(1-2 ± 0-3) X 10^ 
(8-9 ± 1-4) X 10*^ 
(7-2 ± 1-2) X lo" 
(3-7 ± 0-9) X lo" 

(5-7 ± 1-7) X 10-1 



(Hm) log(Af/Mp)) = 9-0-10-0 log(M/M0) = 10-0-11-0 



z = 10-15 

log(M/M0) = 11 



0-12-0 log(M/M0) = 9-0-9-5 log(M/M0 ) = 9-5-10-0 



70 
160 

250 
350 
500 
1100 



(2-8 ± 0-2) X 10" 
(5-1 ± 5-3) X 10" 
(5-0 ± 3-2) X 10" 



(5-5 ± 8-1) : 



(1-4 ± 0-1) 



10' 



(1-3 ± 0-1) X 10" 
(5-8 ± 1-0) X 10-1 
(2-4 ± 1-8) X IQ-^ 



(9-6 ± 0-5) X 1 
(5-4 ± 0-8) X 1 
(2-7 ± 0-5) X 
(4-2 ± 0-2) X 
(3-7 ± 0-2) X 
(2-0 ± 0-2) X 
(1-4 ± 0-3) X 1 



(1-6 ± 0-1) X 10- 
(1-2 ±0-1) X 10*^ 
(6-9 ± 0-8) X lo" 
(8-8 ± 0-4) X lo" 
(8-0 ± 0-4) X 10° 
(4-5 ± 0-3) X 10° 

(4-8 ± 0-5) X 10- 



(2-8 ± 0-2) X 10-^ 
(1-6 ± 0-3) X 10° 
(1-2 ± 0-2) X lol 
(1-4 ± 0-1) X lol 
(1-2 ± 0-1) X lol 
(7-4 ± 0-8) X 10° 

(4-1 ± 1-0) X 10-1 



(urn) log(Af/M0) = 9-0-10-0 log(-i\f/M0) = 10-0-11-0 



z = 1-5-2-0 

log(M/M0) = 11-0-12-0 log(Af/M0) = 9-0-9-5 log(A^/M0 ) = 9-5-10-0 



24 

70 
160 
250 
350 
500 
1100 



(1-8 ± 0-1) X 10-^ 

(2-4 ± 3-9) X 10-1 
(S-6 ± 1-3) X 10-1 
(8.6 ± 1-4) X 10-1 
(6-0 ± 1-3) X 10-1 
(1-6 ± 2-2) X 10-2 



(9-0 ± 0-4) X li 
(2-9 ± 0-9) X 1' 
(1-8 ± 0-6) X : 
(3-1 ± 0-2) X : 
(3-0 ± 0-2) X : 
(1-9 ± 0-2) X : 
(1-5 ± 0-4) X 1' 



(1-8 ± 0-1) X 10-^ 
(4-3 ± 1-3) X 10-1 
(2-9 ± 0-9) X 10° 
(6-5 ± 0-4) X 10° 
(6-4 ± 0-4) X 10° 
(4-5 ± 0-3) X 10° 
(3-7 ± 0-6) X IQ-l 



(3-1 ± 0-3) X 10 

(1-7 ± 0-2) X 10' 

(7-8 ± 1-6) X 10' 

(1-3 ± 0-1) X 10 

(1-3 ± 0-1) X 10 

(8-9 ± 0-7) X 10' 
(1-1 ± 0-1) 



A 

(urn) log(Af/M0) = 9-0-10-0 



log(Af/M0) = 10-0-11-0 
T 



z = 2-0-2-5 

g(A-J/M0) = 11 



0-12-0 
5^ 

0-1 



og(Af/M0) = 9-0-9-5 
(1-3 ± 0-1) X 10-1 
(4-4 ± 1-3) X 10-1 
(2-9 ± 0-9) X 10° 
(5-2 ± 0-3) X 10° 
(6-0 ± 0-4) X 10° 
(4-5 ± 0-3) X 10° 
(4-7 ± 0-7) X 10-1 



log(Af/M0) 



^ 9-5-10-0 



(2-S ± 0-1) X 10-1 
(1-4 ± 0-2) X 10° 
(6-5 ± 1-3) X 10° 
(1-2 ± 0-1) X lol 
(1-3 ± 0-1) X lol 
(1-0 ± 0-1) X loi 



70 
160 
250 
350 
500 
1100 



(1-2 ± 1-6) X 10" 



(7-5 ± 15-0) X 10" 
(2-8 ± 1-6) X 10- 
(3-5 ± 1-7) X 10- 
10" 



(1-9 ± 33-0) 



(5-5 ± 0-4) X 1 
(1-0 ± 1-0) X 
(1-8 ± 0-6) X 
(2-0 ± 0-2) X 
(2-3 ± 0-2) X 
(2-0 ± 0-2) X 
(1-5 ± 0-4) X 



(1-4 ± 0-1) X 10 



A 



(urn) log(Af/MQ) = 9-0-10-0 log(Af/M0 ) = 10-0-11-0 



z = 2-5-3-0 

log(Af/M0) = 11 



0-12-0 log(Af/M0) = 9-0-9-5 log(Af/M0) = 9-5-10-0 



24 
70 
160 

250 
350 
500 
1100 



(9-7 ± 16.0) X 10" 



(2-9 ± 1-5) X 10" 
(4-1 ± 1-7) X 10" 
(3-0 ± 1-7) X 10" 
(4-4 ± 3-7) X 10" 



(2-1 ± 0-2) X 1 
(1-5 ± 1-0) X 1 



(1-4 ± 0-2) X 10° 
(1-6 ± 0-2) X 10° 
(1-0 ± 0-2) X 10° 
(1-6 ± 0-4) X 10-1 



(7-7 ± 0-4) X 10-^ 
(2-3 ± 1-4) X 10-1 
(3-4 ± 0-9) X 10° 
(3-7 ± 0-3) X 10° 
(4-7 ± 0-3) X 10° 
(3-7 ± 0-3) X 10° 
(6-0 ± 0-7) X 10-1 



(2-1 ± 0-2) X 10 
(3-5 ± 2-2) X 10 
(3-1 ± 1-3) X 10' 
(1-1 ± 0-1) X 10 
(1-4 ± 0-1) X 10 
(1-1 ± 0-1) X 10 
(1-9 ± 0-1) X 10' 



(urn) log(A//M0) = 9-0-10-0 log(Af/M0) = 10-0-11-0 1( 



z = 3-0-3-5 

g(A//M0) = 11-0-12-0 



log(Af/M0) = 9-0-9-: 
(6-4 ± 0-5) X 10-2 



log(Af/M0) = 9-5-10-0 



(1-8 ± 0-1) X 10-1 

(1-3 ± 0-3) X 10° 

(1-1 ± 0-2) X lol 

(1-2 ± 0-1) X lol 

(1-5 ± 0-1) X lol 

(1-3 ± 0-1) X lol 

(2-2 ± 0-2) X 10° 



70 
160 
250 
350 
500 
1100 



(3-0 ± 1-5) X 10" 



(3-1 ± 47-0) X 10" 



(8-5 ± 25-0) X 10-2 
(3-1 ± 2-6) X 10-1 



(1-2 ± 0-2) X 10" 

(1-5 ± 0-7) X lo' 
(1-3 ± 0-2) X lo' 
(1-8 ± 0-3) X lo' 
(1-8 ± 0-3) X lo' 
(2-1 ± 0-5) X 10" 



(3-5 ± 1-3) X 10 
(5-1 ± 0-4) X 10' 
(7-1 ± 0-5) X 10' 
(6-6 ± 0-5) X 10' 
(8-7 ± 1-2) X 10" 



n° 



A 



(urn) log(Af/MQ) = 9-0-10-0 log(Af/M0) = 10-0-11-0 



z = 3-5-4-0 

log(-M/M0) = 11-0-12-0 



log(A;/M0) = 9-0-9-5 log(Af/M0) = 9-5-10-0 



70 
160 
250 
350 
500 
1100 



(6-6 ± 0-2) X 10- 
8-5 ± 1000-0) X 10 

(3-3 ± 96-0) X lo' 
(7-0 ± 0-6) X 10° 
(7-0 ± 0-6) X 10° 
(6-2 ± 0-7) X 10° 

(9-9 ± 0-5) X 10- 



(9-9 ± 7-1) X 10" 
(1-3 ± 1-4) X 10" 



(5-7 ± 18-0) X 10-1 
(5-6 ± 5-3) X 10-1 
(7-1 ± 6-3) X 10-1 
(1-0 ± 0-6) X 10° 
(2-3 ± 1-3) X 10-1 



(3-2 ± 1-0) X 10" 
(4-1 ± 4-5) X 10" 



(3-2 ± 1-0) X 10 
(3-8 ± 1-2) X 10' 
(3-3 ± 1-1) X 10' 
(6-4 ± 2-5) X 10" 



(1-1 ± 0-3) X 10 
(3-6 ± 8-7) X 10 
(1-4 ± 5-9) X 10' 
(1-1 ± 0-2) X 10 
(1-6 ± 0-3) X 10 
(1-5 ± 0-2) X 10 
(2-7 ± 0-0) X 10' 



Note- — Stacked Flux 
adequately converge- Fit 



Dens 
X den 



extended bootstrap te 
Figure [7] 
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TABLE 5 
Stacked Flux Densities of Quiescent Galaxies 



Note. — StacI 
adequately conv 



X 




z = 0.0-0.5 






(Hm) 


log(Af/MQ) = 9.0-10.0 


log(*f/M0) = 10.0-11.0 


log(A.J/MQ) = 


11.0-12.0 


24 


(1.5 ± 0.5) X 10-^ 


(2.9 ± 0.4) X 10-^ 


(1.1 ± 0.2) : 


K 10-1 


70 


(9.1 ± 20.0) X 10^^ 


(1.7 ± 1.9) X 10-1 


— 




160 


— 


(1.1 ± 1.1) X lO" 


— 




250 


— 


(5.6 ± 2.8) X 10-1 


(2.4 ± 1.1) 


X lo" 


350 


— 


(2.6 ± 3.0) X 10-1 


(1.3 ± 1.3) 


X lO" 


500 


— 


— 


(7.7 ± 13.0) 


X 10-1 


1100 


— 


— 


— 




\ 




z = 0.5-1.0 






(Hm) 


log(M/MQ) = 9.0-10.0 


lag(M/MQ-j = 10.0-11.0 


lag(M/MQ) = 


11.0-12.0 


24 


— 


(1.4 ± 0.2) X 10--^ 


(2.5 ± 0.7) : 


K 10-^ 


70 


(7.1 ± 11.0) X 10^2 


— 


— 




160 


(9.9 ± 56.0) X 10"^ 


— 


(1.8 ± 1.7) 


X lO" 


250 


(3.8 ± 1.9) X 10"-^ 


— 


(4.7 ± 5.1) : 


< 10-1 


350 


(5.7 ± 2.1) X 10-1 


— 


(2.1 ± 6.0) : 


< 10-1 


500 


(3.6 ± 1.9) X 10-1 


— 


— 




1100 


(4.9 ± 3.6) X 10-^ 


— 


(6.6 ± 11.0) 


X 10-2 


X 




z = 1.0-1.5 






(Hm) 


log(Af/MQ) = 9.0-10.0 


log(*f/M0) = 10.0-11.0 


log(A.J/MQ) = 


11.0-12.0 


24 


— 


(6.2 ± 1.7) X 10--^ 


(2.2 ± 0.6) : 


K 10-2 


70 


— 


— 


(2.6 ± 2.2) : 


< 10-1 


160 


— 


— 


— 




250 


— 


(2.4 ± 1.4) X 10-1 


(1.1 ± 0.4) 


X lo" 


350 


— 


(3.1 ± 1.7) X 10-1 


(1.7 ± 0.5) 


X lO" 


500 


— 


(1.2 ± 1.6) X 10-1 


(1.4 ± 0.5) 


X lo" 


1100 


— 


— 


— 




X 




z = 1.5-2.0 






(Hm) 


log(Af/Mp,) = 9.0-10.0 


log(Af/MQ) = 10.0-11.0 


log(Af/M0) = 


11.0-12.0 


24 


— 


(1.3 ± 1.8) X 10-3 


(1.1 ±0.7) : 


K 10-2 


70 


— 


— 


— 




160 


— 


(2.2 ± 0.7) X 10*^ 


— 




250 


— 


(1.6 ± 1.9) X 10-1 


— 




350 


(2.3 ± 3.4) X 10^1 


(3.9 ± 2.4) X 10-1 


(4.7 ± 7.4) : 


< 10-1 


500 


(3.8 ± 3.6) X 10^1 


(3.7 ± 2.3) X 10-1 


(4.9 ± 7.6) : 


K 10-1 


1100 


— 


— 


(1.7 ± 1.2) : 


K 10-1 


X 




z = 2.0-2.5 






(Hm) 


log(A///Mfr,) = 9.0-10.0 


log(A<f/M0) = 10.0-11.0 


log(Al-/M0) = 


11.0-12.0 


24 


— 


(9.1 ± 3.4) X 10-3 


(6.0 ± 0.7) : 


K 10-2 


70 


(2.8 ± 32.0) X 10-^ 


— 


(3.4 ± 2.8) : 


K 10-1 


160 


(4.8 ± 1.9) X lo" 


(2.2 ± 1100.0) X 10-3 


(3.9 ± 17.0) 


X 10-1 


250 


(9.3 ± 49.0) X 10-^ 


(7.1 ± 3.1) X 10-1 


(4.1 ± 0.6) 


X lO" 


350 


(1.7 ± 0.6) X lo" 


(1.0 ± 0.4) X lo" 


(4.6 ± 0.7) 


X lO" 


500 


(1.9 ± 0.6) X lo" 


(S.2 ± 3.5) X 10-1 


(3.6 ± 0.6) 


X lO" 


1100 


— 


— 


(2.6 ± 1.1) : 


K 10-1 


A 




z = 2.5-3.0 






(Hm) 


log(Af/MQ) = 9.0-10.0 


log(A//MQ) = 10.0-11.0 


log(AJ/MQ) = 


11.0-12.0 


24 


— 


(2.1 ± 0.3) X 10-^ 


(8.5 ± 0.5) : 


.10-2 


70 


— 


— 


(8.7 ± 24.0) 


X 10-2 


160 


— 


(9.2 ± 13.0) X 10-1 


(2.1 ± 1.4) 


X lO" 


250 


— 


(1.1 ± 0.3) X lo" 


(3.8 ± 0.4) 


X lo" 


350 


— 


(1.8 ± 0.4) X lo'^ 


(4.1 ± 0.6) 


X lO" 


500 


— 


(1.8 ± 0.4) X lo'^ 


(3.2 ± 0.5) 


X lo" 


1100 


— 


(2.2 ± 0.8) X 10-1 


(3.6 ± 1.1) : 


< 10-1 


X 




z = 3 0-3 5 






(Hm) 


log(A//MQ) = 9.0-10.0 


log(A//MQ) = 10.0-11.0 


log(Af/MQ) = 


11.0-12.0 


24 


— 


(2.1 ± 0.4) X 10-2 


(8.9 ± 0.7) : 


< 10-2 


70 


(2.7 ± 200.0) X lo" 


— 


— 




160 


(5.2 ± 180.0) X 10*^ 


— 


(5.0 ± 1.8) 


X lO" 


250 


(2.3 ± 1.1) X lo'^ 


(9.9 ± 4.4) X 10-1 


(5.7 ± 0.6) 


X lo" 


350 


(2.5 ± 1.1) X lo" 


(1.6 ± 0.6) X lo" 


(6.2 ± 0.8) 


X lO" 


500 


(4.8 ± 1.3) X lo" 


(1.2 ± 0.5) X lo" 


(4.8 ± 0.8) 


X lofl 


1100 


— 


(2.6 ± 1.1) X 10-1 


(4.6 ± 1.4) : 


K 10-1 


X 




z = 3.5-4.0 






(Hm) 


log(Af/MQ) = 9.0-10.0 


log(Af/M0) = 10.0-11.0 


log(A-J/MQ) = 


11.0-12.0 


24 


— 


(4.4 ± 5.1) X 10-3 


(5.7 ±1.1) ; 


K 10-2 


70 


— 


(1.8 ± 4.6) X 10-1 


— 




160 


— 


(8.2 ± 27.0) X 10-1 


(7.1 ± 3.3) 


X loO 


250 


— 


(3.7 ± 0.8) X lo" 


(3.8 ± 1.1) 


X lo" 


350 


— 


(4.7 ±1.1) X 10*^ 


(5.9 ± 1.4) 


X lO" 


500 


— 


(3.2 ± 1.1) X lo'^ 


(5.9 ± 1.5) 


X lo" 


1100 


— 


(4.7 ± 3.4) X 10-1 


(1.9 ± 0.4) 


X lO" 


aijy. Err 


ora estimated with a extcr 


ided bootstrap technique des 


eribed in § [TI] 


Blanlt space 


onverted 


to ul^, are shown in Figu 


reEl 
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TABLE 6 
Stacked Flux Densities in Bins of Redshift 



Band Type 



500 



z = 0.0-1.0 
(nWm~' 



Stack (7.8 ±0.5) X 10"^ [50%] 
(7.8 ±0.5) X 10"^ [40%] 
.„, , ..N ,, .-^-1 [49%] 



z = 1.0-2.0 
(nWm^^sr^i) 



Stack 
70 CC 

SED 



Stack 
160 CC 

SED 



Stack 
250 CC 

SED 



Stack 
350 CC 

SED 



2.5 ±0.3 
2.5 ±0.3 
2.4 ±0.9 



4.3 ±0.6 
4.3 ±0.6 
4.8 ± 1.8 



2.7 ±0.4 
2.7 ±0.4 
2.9 ± 1.1 



1.3 ±0.3 
1.3 ±0.3 
1.2 ±0.5 



Stack 
CC 
SED 



(4.6 ± 1.0) 

(4.6 ± i.o; 

(3.6 ± 1.5; 



Stack 
1100 CC 
SED 



(1.0 ± 1.4) 
(1.0± 1.4 
(1.4±0.5 



: 10" [66%] 
X 10" [62%] 
X 10" [60%] 



(5.7 ±0.7) X 10^1 [36%] 
(6.2 ± 1.1) X 10"^ [31%] 
(7.3 ±2.7) X 10~^ [37%] 



2 = 2.0-3.0 
(nWm"" 



24 CC (7.8 ± 0.5) X 10"^ [40%] (6.2 ± 1.1) x 10"^ [31%] (3.0 ± 1.4) x 10"^ [15%] 

SED (9.7 ±3.6) X 10 "^ '"""'' .^ o 1 o ^^ ., ,„-l ro^rwi ,n ^ , , n^ .. ,r,-l r,orwi 



10" [49%] 
10" [28%] 
10" [32%] 



(9.9 ±4.3) X 10"^ [26%] 
(9.9 ±6.5) X 10"^ [24%] 
(1.8 ±0.7' " ' -■■ 



(1.8 ±0.3) X 10"^ [11%] 



(3.5 ± 1.3) X 10"^ [18%] 



z = 3.0-4.0 
(nWm^^sr^i) 



(2.5 ±0.7) X 10-2 [1%] 
(2.5 ± 1.0) X 10"^ [12%] 
(3.5 ±0.0) X 10"^ [17%] 



(2.4 ±1.6) X 10"^ [6%] (3.1 ±5.2) x 10"^ [0%] 



(3.1 ±0.9 
(3.3 ± 1.4 
(3.8 ± 1.4 



X 10' 



[39%] 
10° [24%] 
10" [26%] 



(2.9 ±0.7 
(3.2 ± 1.0 
(3.9 ± 1.4 



c 10" [30%] 

c 10" [15%] 

: 10° [14%] 

X 10"^ [23%] 

X 10"^ [9%] 

X 10"^ [7%] 



(2.0 ±0.5 
(2.2 ±0.8 
(2.2 ±0.8 



(8.7 ± 1.6) 



(4.9 ± 12.0) X 10 

f'S S -t- -^ '^'l y 10 



12%] (3.2 ± 150.0) X 10" 



[0%] 



10° [44%] (8.8 ± 3.2) X 10"^ [21%] (9.2 ± 0.0) x 10"^ [22%] 

r. — : Trz : : 1 — : ttt ■ ^ t — : — _ .. 



x 10" [35%] 
x 10° [22%] 
10° [25%] 



(9.1 ±3.6) X 10"^ [10%] 



Tc? 



(1.1 ±2.8) X lO'^ p,„j v„.^^„..y„.„ ^,.,„j 

(1.9 ± 0.7) X 10° [12%] (2.0 ± 0.0) x 10° [13%] 

■ 1 r, .rwi ,„ ^ I „ „N . . ,„-1 r„rwi 



[7%] 

"12%] 



(3.9 ±1.2) X 10-^ [4%] 
(6.2 ±3.1) X 10° [41%] 



[42%] (9.7 ±2.5) X 10-^ [14%] (2.5 ± 0.6) x lO^^ [3%] 



: 10° [29%] (1.8 ± 1.4) x 10° [16%] 

10° [35%] (2.6 ± 0.9) X 10° [23%] 

. _/i r . „^, ttt: — : — : — rr 1 — r. _,w 



(3.3 ± 1.0) X 10° [30%] 
(3.0 ± 0.0) X 10° [27%] 



X 10° [45%] (8.5 ± 1.8) X 10-^ [19%] (2.3 ± 0.5) x lO^^ [5%] 

X 10° [25%] (1.9 ± 1.0) X 10° [21%] (3.1 ± 0.8) x 10° [36%] 

X 10° [26%] (2.0 ± 0.7) X 10° [24%] (3.1 ± 0.0) x 10° [36%] 



.6 ±2.3) 
,2 ± 2.i 



X 10-^ [14%] 
X 10-2 [5%] 
X 10-2 [7%] 



10"^ [43%] (4.9 ±0.6) X 10"^ [24%] 
.„_, ,..„„, ,..,.,., ^po j24%] 



10"^ [20%] 
10-' [17%] 



,7 ±2.3) X 10-2 [38%] 
.8 ±3.4) X 10-2 [15%] 
.8± 1.4) X 10-2 [26%] 



(1.6 ±0.2) X 10"' [8%] 
(1.2 ± 0.4) X 10° [24%] (2.2 ± 0.4) x 10° [45%] 

(9.5 ± 3.2) X 10-'- [20%] (1-8 ± 0-0) x 10° [37%] 



(2.2 ±0.8) X 10" 
(3.7 ±4.5) X 10" 
(5.5 ± 1.8) X 10" 



2 [31%] 



(1.0 ±0.2) X 10-2 [14%] 
[20%] (1.1 ±0.4) X 10-' [58%] 
[30%] (1.3 ± 0.0) X 10-' [73%] 



Note. — Stacked Flux Densities in nijy. Uncorrected flux densities are labeled "Stack", eoniplctcness-correetcd values are labeled "CC", and 
the best-fit SED values are labeled "SED". Errors estimated with a extended bootstrap technique described in ij |3.4l Blank spaces represent where 
the algorithm fails to adequately converge. 



TABLE 7 
Stacked Flux Densities in Bins of Stellar Mass of Star-Forming Galaxies 



Band 


Type 


log(M/MQ) = 9.0-9.5 


log(Af/MQ) = 9.5-10.0 


log(M/MQ 


= 10. 


D-10.5 


log(M/MQ) = 10.5-11.0 


log(A//MQ 


) 


= 11.0-12.0 


(Hm) 




(nWm-2„-l) 


(nWm-2„-l) 


(nWm 


-2ar- 


') 


(„Wm-2sr-l 


) 


(„Wm 


-2„-l) 




Stack 


(9.8 ± 4.6) X 10~2 [e%] 


(3.4 ± 0.4) X 10-1 [21%] 


(5.8 ± 0.3) 


X 10-1 


[37%] 


(3.9 ± 0.1) X 10-1 


[24%] 


(9.5 ± 0.4) 


X 


10-2 [6%] 


24 


CC 


(9 


8 ± 6.7) X 10-2 [5%] 


(4.1 ± 1.5) X 10-1 pjy^] 


(8.9 ± 1.4) 


X 10-1 


[45%] 


(4.0 ± 0.1) X 10-1 


[20%] 


(9.5 ± 0.4) 


X 


10-2 [4%] 




SED 


(1 


3 ± 0.8) X 10-1 [6%] 


(5.4 ± 2.4) X 10-1 [27^^] 


(1.1 ± 0.4) 


X 10" 


[55%] 


(4.7 ± 1.9) X 10-1 


[23%] 


(9.6 ± 4.0) 


X 


10-2 [4%] 




Stack 


(2 


3 ± 3.0) X 10-1 [6%] 


(1.0 ± 0.3) X 10° [27%] 


(1.5 ± 0.2) 


X 10" 


[38%] 


(8.1 ± 0.9) X 10-1 


[21%] 


(1.9 ± 0.2) 


X 


10-1 [6%1 


70 


CC 


(2 


3 ± 4.4) X 10-' [5%] 


(1.0 ± 1.2) X 10° [25%] 


(1.7 ± 1.6) 


X 10" 


[42%] 


(8.2 ± 0.9) X 10-1 


[20%] 


(2.0 ± 0.2) 


X 


10-1 [4%] 




SED 


(3 


5 ± 2.0) X 10-1 [8%] 


(1.4 ± 0.6) X 10° [34%] 


(2.8 ± 1.1) 


X 10" 


[68%] 


(1.1 ± 0.5) X 10° 


27%] 


(2.3 ± 0.9) 


X 


10-1 [6%] 




Stack 


(3 


5 ± 6.2) X 10-1 [4%] 


(1.6 ± 0.6) X 10° [18%] 


(3.4 ± 0.4) 


X 10" 


[38%] 


(2.3 ± 0.2) X 10° 


26%] 


(6.6 ± 0.5) 


X 


10-1 [6%] 


160 


CC 


(3 


5 ± 9.4) X 10-' [2%] 


(1.7 ± 2.9) X 10° [11%] 


(9.2 ± 3.3) 


X 10° 


[62%] 


(2.4 ± 0.2) X 10° 


16%] 


(6.6 ± 0.6) 


X 


10-1 [3%] 




SED 


(6 


9 ± 4.1) X 10-' [4%] 


(2.8 ± 1.2) X 10° [19%] 


(5.9 ± 2.4) 


X 10° 


[39%] 


(2.3 ± 1.0) X 10° 


15%] 


(6.0 ± 2.0) 


X 


10-1 [3%] 




Stack 


(3 


8 ± 4.5) X 10-' [5%] 


(1.4 ± 0.4) X 10° [20%] 


(2.4 ± 0.2) 


X 10° 


[35%] 


(1.8 ± 0.1) X 10° 


26%] 


(6.1 ± 0.2) 


X 


10-1 [7%] 


250 


CC 


(<1 


1 ± 6.6) X 10-' [3%] 


(1.8 ± 1.5) X 10° [16%] 


(6.0 ± 1.4) 


X 10° 


[54%] 


(1.9 ± 0.1) X 10° 


17%] 


(6.1 ± 0.2) 


X 


10-1 [4%] 




SED 


(4 


1 ± 2.4) X 10-1 [3%] 


(2.3 ± 0.9) X 10° [20%] 


(6.2 ± 2.8) 


X 10" 


[56%] 


(2.4 ± 1.0) X 10° 


22%] 


(6.2 ± 2.6) 


X 


10-1 [5%] 




Stack 


(1 


8 ± 3.3) X 10-1 [4%] 


(8.2 ± 3.0) X 10-1 [18%] 


(1.5 ± 0.2) 


X 10" 


[33%] 


(1.3 ± 0.1) X 10° 


28%] 


(3.9 ± 0.2) 


X 


10-1 [8%] 


350 


CC 


(1 


8 ± 4.9) X 10-1 [2%] 


(1.5 ± 1.1) X 10° [17%] 


(4.8 ± 1.0) 


X 10" 


[56%] 


(1.4 ± 0.1) X 10° 


16%] 


(4.0 ± 0.2) 


X 


10-1 [4%] 




SED 


(1 


9 ± 1.0) X 10-1 [2%] 


(1.3 ± 0.5) X 10° [15%] 


(4.8 ± 2.6) 


X 10° 


[56%] 


(1.5 ± 0.6) X 10° 


17%] 


(4.6 ± 1.9) 


X 


10-1 [6%] 




Stack 


(7.2 ± 10.0) X 10-2 [3%] 


(3.2 ± 0.9) X 10-1 [16%] 


(6.4 ± 0.6) 


xlO-1 


[32%] 


(6.0 ± 0.3) X 10-1 


[30%] 


(2.1 ± 0.1) 


X 


10-1 jj^fljj] 


500 


CC 


(7.2 ± 16.0) X 10-2 [i%] 


(8.3 ± 3.8) X 10-1 [17%] 


(2.8 ± 0.4) 


X 10° 


[59%] 


(6.7 ± 0.3) X 10-1 


[14%] 


(2.1 ± 0.1) 


X 


10-1 [4%] 




SED 


(6.6 ± 3.0) X 10-2 [1%] 


(5.5 ± 2.0) X 10-1 [11%] 


(2.3 ± 1.5) 


X 10" 


[49%] 


(6.1 ± 2.4) X 10-1 


[12%] 


(2.1 ± 0.9) 


X 


10-1 [4%] 




Stack 


(7.6 ± 1500.0) X 10-^ [0%] 


(5.0 ± 13.0) X 10-3 [7%] 


(2.2 ± 0.8) 


XlO-2 


[31%] 


(2.6 ± 0.4) X 10-2 


[36%] 


(1.2 ±0.1) 


X 


10-2 JJ75J] 


1100 


CC 


(7.6 ± 2200.0) X 10~^ [0%] 


(6.9 ± 49.0) X lO-'^ [3%] 


(1.3 ± 0.6) 


X 10-1 


[70%] 


(2.8 ± 0.4) X 10-2 


[15%] 


(1.3 ± 0.1) 


X 


10-2 [6%] 




SED 


(3.5 ± 1.1) X 10-^ [1%] 


(3.2 ± 1.2) X 10-2 [17%] 


(1.6 ± 1.1) 


X 10-1 


[S4%] 


(3.1 ± 1.2) X 10-2 


[17%] 


(1.3 ± 0.5) 


X 


10-2 [7%] 



Note. — Stacked 
labeled "SED". Bla 



^ labeled "Sta 



nk spaces represent ' 



fails to adequately 1 



are labeled "CC 
ed bootstrap tec 



d the 
hnique dc; 



:st-fit SEDvah 
■ibed in g f3T4l 
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TABLE 8 

Stacked Flux Densities in Bins of Stellar Mass of Quiescent Galaxies 



Band 


Type 


log(M/MQ) = 9.0-10.0 log(M/M0) = 10.0 


-11.0 log(Af/M£ 


) = 11.0-12.0 


(yxm) 




(nWm-2sr-i) 


(nW m 2 SI- 1 




(nWm-2sr-i) 




Stack 


(1.7 ±4.0) X 10-^ [0%] 


(3.5 ±1.5) X 10-2 


[2%] 


1.8 ±0.2 


X 10-2 [1%] 


24 


CC 


(1.9 ±4.3) X 10-=* [0%] 


(3.8 ±1.6) X 10-2 


[1%] 


1.8 ±0.2 


X 10-2 [0%] 




SED 


(6.6 ± 0.0) X 10-3 [0%] 


(4.6 ±1.8) X 10-2 


[2%] 


2.0 ±0.8 


X 10-2 [1%] 




Stack 


(1.7 ±3.6) X 10-2 [0%] ( 


1.6 ±11.0) X 10-2 


[0%] 


2.1 ±2.4 


X 10-2 [0%] 


70 


CC 


(1.9 ±3.8) X 10-2 [0%] ( 


1.8 ±12.0) X 10-2 


[0%] 


2.1 ±2.4 


X 10-2 [0%] 




SED 


(1.7 ± 0.0) X 10-2 [o%] 


(1.2 ±0.4) X 10-1 


[2%] 


5.2±2.1 


X 10-2 [1%] 




Stack 


(1.6 ±8.1) X 10-2 [0%] 


(3.5 ±2.5) X 10-1 


[3%] 


1.2 ±0.6 


X 10-1 [1%] 


160 


CC 


(1.7 ±8.7) X 10-2 [0%] 


(3.7 ±2.7) X 10-1 


[2%] 


1.2 ±0.6 


X 10-1 [0%] 




SED 


(3.1 ±0.0) X 10-2 [0%] 


(2.3 ±0.9) X 10-1 


[1%] 


1.1 ±0.4 


X 10-1 [0%] 




Stack 


(1.7 ±4.0) X 10-2 [0%] 


(1.3 ±1.5) X 10-1 


[1%] 


1.0 ±0.2 


X 10-1 [1%] 


250 


CC 


(1.9 ±4.3) X 10-2 [0%] 


(1.4 ±1.6) X 10-1 


[1%] 


1.0 ±0.2 


X 10-1 [0%] 




SED 


(3.6 ± 0.0) X 10-2 [0%] 


(2.1 ±0.7) X 10-1 


[1%] 


1.3 ±0.6 


X 10-1 [1%] 




Stack 


(2.1 ±3.0) X 10-2 [0%] 


(1.3 ±1.1) X 10-1 


[3%] 


,8.8 ±1.6 


X 10-2 [1%] 


350 


CC 


(2.3 ±3.2) X 10-2 [0%] 


(1.4 ±1.2) X 10-1 


[1%] 


8.8 ±1.6 


X 10-2 [1%] 




SED 


(3.0 ± 0.0) X 10-2 [0%] 


(1.5 ±0.6) X 10-1 


[1%] 


1.0 ±0.4 


X 10-1 [1%] 




Stack 


(1.1 ± 1.1) X 10-2 [0%] 


(7.2 ±3.6) X 10-2 


[3%] 


5.1 ±0.7 


X 10-2 [2%] 


500 


CC 


(1.2 ± 1.2) X 10-2 [0%] 


(7.8 ±3.9) X 10-2 


[1%] 


5.2 ±0.7 


X 10-2 [1%] 




SED 


(1.6 ±0.0) xlO-2 [0%] 


(8.0 ±2.8) X 10-2 


[1%] 


4.8 ±1.7 


X 10-2 [0%] 




Stack 


(4.9 ± 14.0) x 10-4 [0%] 


(2.3 ±5.0) X 10-3 


[3%] 


2.5 ±0.8 


X 10-3 [3%] 


1100 


CC 


(5.2 ± 15.0) X lO--* [0%] 


(2.4 ±5.4) X 10-3 


[1%] 


2.6 ±0.8 


X 10-3 [1%] 




SED 


(2.7 ±0.0) X 10-3 [1%] 


(1.1 ±0.3) X 10-2 


[6%] 


3.3 ±1.2 


X 10-3 [1%] 



Note. — Stacked flux densities in mjy. Uncorrected flux densities are labeled "Stack", completeness-corrected values are labeled "CC", 
and the best-fit SED values are labeled "SED" . Bla nk s paces represent where the algorithm fails to adequately converge. Errors estimated 
with a extended bootstrap technique described in § 13.41 



